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ABSTRACT 

We review and complete the kinetic theory of spatially inhomogeneous stellar systems when collective effects (dressing 
of the stars by their polarization cloud) are neglected. We start from the BBGKY hierarchy issued from the Liouville 
equation and consider an expansion in powers of 1/N in a proper thermodynamic limit. For TV — > +oo, we obtain 
the Vlasov equation describing the evolution of collisionless stellar systems like elliptical galaxies. At the order 1/N, 
we obtain a kinetic equation describing the evolution of collisional stellar systems like globular clusters. This equation 
coincides with the generalized Landau equation derived by Kandrup (1981) from a more abstract projection operator 
formalism. This equation does not suffer logarithmic divergences at large scales since spatial inhomogeneity is explicitly 
taken into account. Making a local approximation, and introducing an upper cut-off at the Jeans length, it reduces 
to the Vlasov-Landau equation which is the standard kinetic equation of stellar systems. Our approach provides a 
simple and pedagogical derivation of these important equations from the BBGKY hierarchy which is more rigorous 
for systems with long-range interactions than the two-body encounters theory. Making an adiabatic approximation, we 
write the generalized Landau equation in angle-action variables and obtain a Landau-type kinetic equation that is valid 
for fully inhomogeneous stellar systems and is free of divergences at large scales. This equation is less general than the 
Lenard-Balescu-type kinetic equation recently derived by Heyvaerts (2010) since it neglects collective effects, but it is 
substantially simpler and could be useful as a first step. We discuss the evolution of the system as a whole and the 
relaxation of a test star in a bath of field stars. We derive the corresponding Fokker-Planck equation in angle-action 
variables and provide expressions for the diffusion coefficient and friction force. 
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1. Introduction 

In its simplest description, a stellar system can be viewed 
as a collection of N classical point mass stars in Newtonian 
gravitational interaction (Binney and Tremaine 2008). As 
understood early by Henon (1964), self-gravitating systems 
experience two successive types of relaxation: A rapid "col- 
lisionless" relaxation towards a quasi stationary state (QSS) 
that is a virialized state in mechanical equilibrium but not 
in thermodynamical equilibrium, followed by a slow "colli- 
sional" relaxation. One might think that, due to the devel- 
opment of stellar encounters, the system will reach, for suffi- 
ciently large times, a statistical equilibrium state described 
by the Maxwcll-Boltzmann distribution. However, it is well- 
known that unbounded stellar systems cannot be in strict 
statistical equilibrium^ due to the permanent escape of high 
energy stars (evaporation) and the gravothermal catastro- 
phe (core collapse). Therefore, the statistical mechanics of 
stellar systems is essentially an out-of-cquilibrium problem 
which must be approached through kinetic theories. 

The first kinetic equation was written by Jeans (1915). 
Neglecting encounters between stars, he described the dy- 
namical evolution of stellar systems by the collisionless 
Boltzmann equation coupled to the Poisson equation. This 
purely mean field description applies to large groups of stars 



such as elliptical galaxies whose ages are much less than the 
collisional relaxation time. A similar equation was intro- 
duced by Vlasov (1938) in plasma physics to describe the 
collisionless evolution of a system of electric charges inter- 
acting by the Coulomb force. The collisionless Boltzmann 
equation coupled sclf-consistcntly to the Poisson equation 
is oftentimes called the Vlasov equation^, or the Vlasov- 
Poisson system. 

The concept of collisionless relaxation was first un- 
derstood by Henon (1964) and King (1966). Lynden-Bell 
(1967) developed a statistical theory of this process and 
coined the term "violent relaxation". In the collisionless 
regime, the system is described by the Vlasov-Poisson sys- 
tem. Starting from an unsteady or unstable initial con- 
dition, the Vlasov-Poisson system develops a complicated 
mixing process in phase space. The Vlasov equation, be- 
ing time-reversible, never achieves a steady state but de- 
velops filaments at smaller and smaller scales. However, 
the coarse-grained distribution function usually achieves a 
steady state on a few dynamical times. Lynden-Bell (1967) 
tried to predict the QSS resulting from violent relaxation 
by developing a statistical mechanics of the Vlasov equa- 
tion. He derived a distribution function formally equiva- 
lent to the Fermi-Dirac distribution (or to a superposition 
of Fermi-Dirac distributions). However, when coupled to 



For reviews about the statistical mechanics of self- 
gravitating systems see, e.g., Padmanabhan (1990), Katz (2003), 
and Chavanis (2006). 



2 See Henon (1982) for a discussion about the name that one 
should give to that equation. 
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the Poisson equation, these distributions have an infinite 
mass. Therefore, the Vlasov-Poisson system has no statisti- 
cal equilibrium state (in the sense of Lyndcn-Bcll) . This 
is a clear evidence that violent relaxation is incomplete 
(Lynden-Bell 1967). Incomplete relaxation is due to ineffi- 
cient mixing and non-ergodicity. In general, the fluctuations 
of the gravitational potential <5$(r,i) that drive the colli- 
sionlcss relaxation last only for a few dynamical times and 
die out before the system has mixed efficiently (Tremainc 
et al. 1986). Understanding the origin of incomplete relax- 
ation, and developing models of incomplete violent relax- 
ation to predict the structure of galaxies, is a very difficult 
problem (Arad & Johansson 2005). Some models of incom- 
plete violent relaxation have been proposed based on differ- 
ent physical arguments (Bcrtin and Stiavelli 1984, Stiavelli 
and Bertin 1987, Hjorth and Madsen 1991, Chavanis 1998, 
Levin et al. 2008). 

On longer timescales, stellar encounters (sometimes re- 
ferred to as "collisions" by an abuse of language) must 
be taken into account. This description is particularly im- 
portant for small groups of stars such as globular clusters 
whose ages are of the order of the collisional relaxation time. 
Chandrasekhar (1942, 1943a, 1943b) developed a kinetic the- 
ory of stellar systems in order to determine the timescale 
of collisional relaxation and the rate of escape of stars from 
globular clusters^. To simplify the kinetic theory, he con- 
sidered an infinite homogeneous system of stars. He started 
from the general Fokker-Planck equation and determined 
the diffusion coefficient and the friction force (first and sec- 
ond moments of the velocity increments) by considering 
the mean effect of a succession of two-body encounters^ 
However, his approach leads to a logarithmic divergence 
at large scales that is more difficult to remove in stellar 
dynamics than in plasma physics because of the absence of 
Debye shielding for the gravitational forcifl Chandrasekhar 
and von Neumann (1942) developed a completely stochas- 
tic formalism of gravitational fluctuations and showed that 
the fluctuations of the gravitational force are given by the 
Holtzmark distribution (a particular Levy law) in which 
the nearest neighbor plays a prevalent role. From these re- 
sults, they argued that the logarithmic divergence has to be 
cut-off at the interparticle distance (see also Jeans 1929). 
However, since the interparticle distance is smaller than 

3 Early estimates of the relaxation time of stellar systems were 
made by Schwarzschild (1924), Rosseland (1928), Jeans (1929), 
Smart (1938), and Spitzer (1940). On the other hand, the evap- 
oration time was first estimated by Ambartsumian (1938) and 
Spitzer (1940). 

4 Later, Cohen et al. (1950), Gasiorowicz et al. (1956), and 
Rosenbluth et al. (1957) proposed a simplified derivation of the 
coefficients of diffusion and friction. 

5 A few years earlier, Landau (1936) had developed a kinetic 
theory of Coulombian plasmas based on two-body encounters. 
Starting from the Boltzmann (1872) equation, and making a 
weak deflection approximation, he derived a kinetic equation 
for the collisional evolution of neutral plasmas. His approach 
leads to a divergence at large scales that he removed heuristi- 
cally by introducing a cut-off at the Debye length (Debye and 
Hiickel 1923) which is the size over which the electric field pro- 
duced by a charge is screened by the cloud of opposite charges. 
Later, Lenard (1960) and Balescu (1960) developed a more pre- 
cise kinetic theory taking collective effects into account. They 
derived a more elaborate kinetic equation, free of divergence at 
large scales, in which the Debye length appears naturally. This 
justifies the heuristic procedure of Landau. 



the Debye length, the same arguments should also apply in 
plasma physics, which is not the case. Therefore, the conclu- 
sions of Chandrasekhar and von Neumann are usually taken 
with circumspection. In particular, Cohen et al. (1950) ar- 
gue that the logarithmic divergence should be cut-off at the 
Jeans length which gives an estimate of the system's size. 
Indeed, while in neutral plasmas the effective interaction 
distance is limited to the Debye length, in a self-gravitating 
system, the distance between interacting particles is only 
limited by the system's size. Therefore, the Jeans length is 
the gravitational analogue of the Debye length. These ki- 
netic theories lead to a collisional relaxation time scaling as 
t r ~ (N/ \nN)tD, where tu is the dynamical time and N 
the number of stars in the system. Chandrasekhar (1949) 
also developed a Brownian theory of stellar dynamics and 
showed that, from a qualitative point of view, the results of 
kinetic theory can be understood very simply in that frame- 
work. In particular, he showed that a dynamical friction is 
necessary to maintain the Maxwcll-Boltzmann distribution 
of statistical equilibrium and that the coefficients of friction 
and diffusion are related to each other by an Einstein rela- 
tion which is a manifestation of the fluctuation-dissipation 
theorem. This relation is confirmed by his more precise ki- 
netic theory based on two-body encounters. It is impor- 
tant to emphasize, however, that Chandrasekhar did not 
derive the kinetic equation for the evolution of the system 
as a whole. Indeed, he considered the Brownian motion of a 
test star in a fixed distribution of field stars (bath) and de- 
rived the corresponding Fokker-Planck equation. This equa- 
tion has been used by Chandrasekhar (1943b), Spitzer and 
Harm (1958), Michie (1963), King (1965), and more re- 
cently Lemon and Chavanis (2010) to study the evaporation 
of stars from globular clusters in a simple setting. 

King (1960) noted that, if we were to describe the dy- 
namical evolution of the cluster as a whole, the distribution 
of the field stars must evolve in time in a self-consistent 
manner so that the kinetic equation must be an integrodif- 
ferential equation. The kinetic equation obtained by King, 
from the results of Rosenbluth et al. (1957), is equivalent to 
the Landau equation, although written in a different form. 
It is interesting to note, for historical reasons, that none 
of the previous authors seemed to be aware of the work 
of Landau (1936) in plasma physics. There is, however, an 
important difference between stellar dynamics and plasma 
physics. Neutral plasmas arc spatially homogeneous due 
to Debye shielding. By contrast, stellar systems are spa- 
tially inhomogeneous. The above-mentioned kinetic theo- 
ries developed for an infinite homogeneous system can be 
applied to an inhomogeneous system only if we make a lo- 
cal approximation. In that case, the collision term is calcu- 
lated as if the system were spatially homogeneous or as 
if the collisions could be treated as local. Then, the ef- 
fect of spatial inhomogeneity is only retained in the ad- 
vection (Vlasov) term which describes the evolution of the 
system due to mean-field effects. This leads to the Vlasov- 
Landau equation which is the standard kinetic equation 
of stellar dynamics. To our knowledge, this equation has 
been first written (in a different form), and studied, by 
Hcnon (1961). Henon also exploited the timescale separa- 
tion between the dynamical time and the relaxation 
time tji 3> tjj to derive a simplified kinetic equation for 
/(e, t), where e = v 2 /2 + <&(r, t) is the individual energy of 
a star by unit of mass, called the orbit-averaged Fokker- 
Planck equation. In this approach, the distribution func- 
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tion /(r,v,i), averaged over a short timescale, is a steady 
state of the Vlasov equation of the form /(e, t) which slowly 
evolves in time, on a long timescale, due to the develop- 
ment of "collisions" (i.e. correlations caused by finite N 
effects or graininess). Henon used this equation to obtain a 
more relevant value for the rate of evaporation from globu- 
lar clusters, valid for inhomogeneous systems. Cohn (1980) 
numerically solved the orbit-averaged Fokkcr-Planck equa- 
tion to describe the collisional evolution of star clusters. His 
treatment accounts both for the escape of high energy stars 
put forward by Spitzer (1940) and for the phenomenon of 
core collapse resulting from the gravothermal catastrophe 
discovered by Antonov (1962) and Lynden-Bell and Wood 
(1968) on the basis of thermodynamics and statistical me- 
chanics. The local approximation, which is a crucial step in 
the kinetic theory, is supported by the stochastic approach 
of Chandrasekhar and von Neumann (1942) showing the 
preponderance of the nearest neighbor. However, this re- 
mains a simplifying assumption which is not easily control- 
lable. In particular, as we have already indicated, the local 
approximation leads to a logarithmic divergence at large 
scales that is difficult to remove. This divergence would 
not have occurred if full account of spatial inhomogeneity 
had been given since the start. 

The effect of spatial inhomogeneity was investigated 
by Severne and Haggcrty (1976), Parisot and Severne 
(1979), Kandrup (1981), and Chavanis (2008). In partic- 
ular, Kandrup (1981) derived a generalized Landau equa- 
tion from the Liouville equation by using projection op- 
erator technics. This generalized Landau equation is in- 
teresting because it takes into account effects of spa- 
tial inhomogeneity which were neglected in previous ap- 
proaches. Since the finite extension of the system is prop- 
erly accounted for, there is no divergence at large scalcfjj. 
Furthermore, this approach clearly show which approxima- 
tions arc needed in order to recover the traditional Landau 
equation. Unfortunately, the generalized Landau equation 
remains extremely complicated for practical applications. 

In addition, this equation is still approximate as it ne- 
glects collective effects and considers binary collisions be- 
tween "naked" particles. As in any weakly coupled system, 
the particles engaged in collisions are "dressed" by the po- 
larization clouds caused by their own influence on other 
particles. Collisions between dressed particles have quan- 
titatively different outcomes than collisions between naked 
ones. In the case of plasmas, collective effects are responsi- 
ble for Debye shielding and they are accounted for in the 
Lenard-Balescu equation. They allow to eliminate the loga- 
rithmic divergence that occurs at large scales in the Landau 
equation. For self-gravitating systems, they lead to "anti- 
shiclding" and are more difficult to analyze because the 
system is spatially inhomogeneoufl If we consider a finite 



There remains a logarithmic divergence at small scales due 
to the neglect of strong collisions. This divergence can be cured 
heuristically by introducing a cut-off at the Landau length which 
corresponds to the impact parameter leading to a deflection at 
90° (Landau 1936). 

7 In a plasma, since the Coulomb force between electrons is re- 
pulsive, each particle of the plasma tends to attract to it particles 
of opposite charge and to repel particles of like charge, thereby 
creating a kind of "cloud" of opposite charge which screens the 
interaction at the scale of the Debye length. In the gravitational 
case, since the Newton force is attractive, the situation is con- 
siderably different. The test star draws neighboring stars into 



homogeneous system, and take collective effects into ac- 
count, one finds a severe divergence of the diffusion coeffi- 
cient when the size of the domain reaches the Jeans scale 
(Weinberg 1993). This divergence, which is related to the 
Jeans instability, does not occur in a stable spatially inho- 
mogeneous stellar system. Some authors like Thorne (1968), 
Miller (1968), Gilbert (1968,1970), and Lerche (1971) at- 
tempted to take collective effects and spatial inhomogeneity 
into account. However, they obtained very complicated ki- 
netic equations that have not found application until now. 
They managed, however, to show that collective effects arc 
equivalent to increasing the effective mass of the stars, 
hence diminishing the relaxation time. Since, on the other 
hand, the effect of spatial inhomogeneity is to increase the 
relaxation time (Parisot and Severne 1979), the two effects 
act in opposite direction and may balance each other. 

Recently, Heyvaerts (2010) derived from the BBGKY 
hierarchy issued from the Liouville equation a kinetic equa- 
tion in angle-action variables that takes both spatial inho- 
mogencities and collective effects into account. To calculate 
the collective response, he used Fourier-Laplace transforms 
and introduced a bi-orthogonal basis of pairs of density- 
potential functions (Kalnajs 1971a). The kinetic equation 
derived by Heyvaerts is the counterpart for spatially inho- 
mogeneous self-gravitating systems of the Lenard-Balescu 
equation for plasmas. Following his work, we showed that 
this equation could be obtained cquivalcntly from the 
Klimontovich equation by making the so-called quasilinear 
approximation (Chavanis 2012). We also developed a test 
particle approach and derived the corresponding Fokkcr- 
Planck equation in angle-action variables, taking collective 
effects into account. This provides general expressions of 
the diffusion coefficient and friction force for spatially in- 
homogeneous stellar systems. 

In a sense, these equations solve the problem of the ki- 
netic theory of stellar systems since they take into account 
both spatial inhomogeneity and collective effects. However, 
their drawback is that they are extremely complicated to 
solve (in addition of being complicated to derive). In an 
attempt to reduce the complexity of the problem, we shall 
derive in this paper a kinetic equation that is valid for spa- 
tially inhomogeneous stellar systems but that neglects col- 
lective effects. Collective effects may be less crucial in stellar 
dynamics than in plasma physics. In plasma physics, they 
must be taken into account in order to remove the diver- 
gence at large scales that appears in the Landau equation. 
In the case of stellar systems, this divergence is removed 
by the spatial inhomogeneity of the system, not by col- 
lective effects. Actually, previous kinetic equations based 
on the local approximation ignore collective effects and al- 
ready give satisfactory results. We shall therefore ignore 
collective effects and derive a kinetic equation (in position- 
velocity and angle-action variables) that is the counter- 
part for spatially inhomogeneous self-gravitating systems of 
the Landau equation for plasmas. Our approach has three 
main interests. First, the derivation of this Landau-type 
kinetic equation is considerably simpler than the deriva- 
tion of the Lenard-Balescu-type kinetic equation, and it can 
be done in the physical space without having to introduce 



its vicinity and these add their gravitational force to that of the 
test star itself. The "bare" gravitational force of the test star is 
thus augmented rather than shielded. The polarization acts to 
increase the effective gravitational mass of a test star. 



3 



Pierre-Henri Chavanis: Kinetic theory of spatially inhomogeneous stellar systems without collective effects 



Laplace- Fourier transforms nor bi-orthogonal basis of pairs 
of density-potential functions. This offers a more physical 
derivation of kinetic equations of stellar systems that can 
be of interest for astrophysicists. Secondly, our approach is 
sufficient to remove the large-scale divergence that occurs 
in kinetic theories based on the local approximation. It rep- 
resents therefore a conceptual progress in the kinetic theory 
of stellar systems. Finally, this equation is simpler than the 
Lenard-Balcscu-type kinetic equation derived by Heyvaerts 
(2010), and it could be useful in a first step before consid- 
ering more complicated effects. Its drawback is to ignore 
collective effects but this may not be crucial as we have ex- 
plained (as suggested by Weinberg's work, collective effects 
become important only when the system is close to insta- 
bility). This Landau- type equation was previously derived 
for systems with arbitrary long-range interaction^ in vari- 
ous dimensions of space (see Chavanis 2007,2008,2010) but 
we think that it is important to discuss these results in the 
specific case of self-gravitating systems with complements 
and amplification. 

The paper is organized as follows. In Section [21 we study 
the dynamical evolution of a spatially inhomogeneous stel- 
lar system as a whole. Starting from the BBGKY hierar- 
chy issued from the Liouville equation, and neglecting col- 
lective effects, we derive a general kinetic equation valid 
at the order 1/N in a proper thermodynamic limit. For 
N —> +oo, it reduces to the Vlasov equation. At the or- 
der 1/N we recover the generalized Landau equation de- 
rived by Kandrup (1981) from a more abstract projection 
operator formalism. This equation is free of divergence at 
large scales since spatial inhomogeneity has been properly 
accounted for. Making a local approximation and introduc- 
ing an upper cut-off at the Jeans length, we recover the 
standard Vlasov-Landau equation which is usually derived 
from a kinetic theory based on two-body encounters. Our 
approach provides an alternative derivation of this funda- 
mental equation from the more rigorous Liouville equation. 
It has therefore some pedagogical interest. In Section [3l we 
study the relaxation of a test star in a steady distribution 
of field stars. We derive the corresponding Fokker-Planck 
equation and determine the expressions of the diffusion and 
friction coefficients. We emphasize the difference between 
the friction by polarization and the total friction (this dif- 
ference may have been overlooked in previous works). For a 
thermal bath, we derive the Einstein relation between the 
diffusion and friction coefficients and obtain the explicit ex- 
pression of the diffusion tensor. This returns the standard 
results obtained from the two-body encounters theory but, 
again, our presentation is different and offers an alterna- 
tive derivation of these important results. For that reason, 
we give a short review of the basic formulae. In Section 21 
we derive a Landau-type kinetic equation written in angle- 
action variables and discuss its main properties. This equa- 
tion, which does not make the local approximation, applies 
to fully inhomogeneous stellar systems and is free of diver- 
gence at large scales. We also develop a test particle ap- 
proach and derive the corresponding Fokker-Planck equa- 
tion in angle-action variables. Explicit expressions are given 
for the diffusion tensor and friction force, and they are com- 
pared with previous expressions obtained in the literature. 



8 For a review on the dynamics and thermodynamics of sys- 
tems with long-range interactions, see Campa et al. (2009). 



2. Evolution of the system as a whole 

2.1. The BBGKY hierarchy 

We consider an isolated system of N stars with identical 
mass m in Newtonian interaction. Their dynamics is fully 
described by the Hamilton equations 



dvi dH dvi 
m = — — , m 

dt 



dt dvi ' 



dH 
dr., ' 



(1) 



i<j 



This Hamiltonian system conserves the energy E = H, 
the mass M = Nm, and the angular momentum L = 
Y2i mr i x v i- As recalled in the Introduction, stellar sys- 
tems cannot reach a statistical equilibrium state in a strict 
sense. In order to understand their evolution, it is necessary 
to develop a kinetic theory. 

We introduce the TV-body distribution function 
Pjv (ri, Vi, ...,rjv, v;v,t) giving the probability density of 
finding, at time t, the first star with position ri and veloc- 
ity Vi, the second star with position r2 and velocity V2 etc. 
Basically, the evolution of the iV-body distribution function 
is governed by the Liouville equation 



dP, 



N 



dt 



where 



N 

E 



dvi 



dP, 



N 



dP> 



N 



dr t 



dvi 



= 0, 



(2) 



(3) 



is the gravitational force by unit of mass experienced by the 
i-th star due to its interaction with the other stars. Here, 
$d(r) denotes the exact gravitational potential produced by 
the discrete distribution of stars and F(j — s- i) denotes the 
exact force by unit of mass created by the j'-th star on the 
i-th star. The Liouville equation ([2]), which is equivalent to 
the Hamilton equations ([1} , contains too much information 
to be exploitable. In practice, we are only interested in the 
evolution of the one-body distribution Pi(r,v,i). 

From the Liouville equation @ we can construct the 
complete BBGKY hierarchy for the reduced distribution 
functions 

Pj(xi, -,Xj,i) = J P N (x 1 ,...,x N ,t)dx j+1 ...dx N , (4) 

where the notation x stands for (r, v). The generic term of 
this hierarchy reads 



dPi 



j 



dP, 



dP n 



-^ + E^+E E F ( fc -*)-^ 

i—l i—l k—l,k^i 

+(n - i) e y p y + 1 o • = °- 



(5) 



This hierarchy of equations is not closed since the equation 
for the one-body distribution Pi (xi , t) involves the two- 
body distribution P2(xi,X2,i), the equation for the two- 
body distribution P-i (xi , X2 , t) involves the three-body dis- 
tribution Ps(xx, X2, X3, t), and so on. 
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It is convenient to introduce a cluster representation of 
the distribution functions. Specifically, we can express the 
reduced distribution Pj(xi, Xj, t) in terms of products 
of distribution functions P,<<j(xi, ...,Xj/,i) of lower order 
plus a correlation function Pj (xi , Xj, t) [see, e.g., Eqs. © 
and © below] . Considering the scaling of the terms in each 
equation of the BBGKY hierarchy, we can see that there 
exist solutions of the whole BBGKY hierarchy such that the 
correlation functions P' scale as 1/N^~ l in the proper ther- 
modynamic limit N — > +00 defined in Appendix [21 This 
implicitly assumes that the initial condition has no corre- 
lation, or that the initial correlations respect this scaling. 
If this scaling is satisfied, we can consider an expansion 
of the BBGKY hierarchy in terms of the small parameter 
1/A. This is similar to the expansion of the BBGKY hi- 
erarchy is plasma physics in terms of the small parameter 
1/A where A> 1 represents the number of charges in the 
Debye sphere (Balescu 2000). However, in plasma physics, 
the system is spatially homogeneous (due to Debye shield- 
ing which restricts the range of interaction) while, for stellar 
systems, spatial inhomogencity must be taken into account. 
This brings additional terms in the BBGKY hierarchy that 
are absent in plasma physics. 



2.2. The truncation of the BBGKY hierarchy at the order 

1/N 

The first two equations of the BBGKY hierarchy are 



^(l) + Vl .— (1) 



+ (N-l) /F(2-H). |^(1,2) dx a =0, (6) 



^(l,2) + v 1 .^(l,2 )+ F(2^1). ^(1,2) 

2 Ot OTi OVi 

/BP 
F(3 1) • -^(1, 2, 3) dx 3 + (1 O 2) = 0. (7) 

We decompose the two- and three-body distributions in the 
form 



P 2 (x x , X 2 ) = Pr (xi )Pi (x 2 ) + P' 2 ( Xl , x 2 ) 



(8) 



P 3 (xi,x 2 ,x 3 ) = Pi(xi)P 1 (x 2 )P 1 (x 3 ) + P 2 '(x 1 ,x 2 )P 1 (x 3 ) 

+P^ ( Xl , X 3 )P X (X 2 ) + P' 2 (X 2 , X 3 )Pl (Xx ) + P^ (Xi , X 2 , X 3 ), 

(9) 

where P'(xi, Xj, t) is the correlation function of order 
j. Substituting Eqs. © and © in Eqs. © and ©, and 
simplifying some terms, we obtain 



-^(1) + vr ■ — (1) 



+ (N-l) 



-(N-l) 



Jf(2^ l)Pi(2)dx 2 



9vi 



yF(2^1)P^(l,2)dx 2 , (10) 



If there are non-trivial correlations in the initial state, like 
primordial binary stars, the kinetic theory will be different from 
the one developed in the sequel. 



J^M + vi-g<l.*> + F(^l).§3 M 



+(N-2) Jf(3^ l)Pi(3)dx 3 
F(2->-l)-yP(3->-l)Pi(3)d3C3 



gP| 
9vi 



1,2) 



gPi 
9vi 



(1)A(2) 



+(A-2) 



/" F(3^ l)P 2 '(2,3)dx 3 



gPi 
9vi 



(1) 



d f 

-— I F(3^1)P^(l,3)P 1 (2)dx 3 

+(N - 2)^- ■ J F(3 -> 1)^(1, 2, 3) dx 3 + (1 o 2) = 0. 

(11) 

Equations (| 10[) and (fTT|) arc exact for all N but they are 
not closed. As explained previously, we shall close these 
equations at the order 1/A in the thermodynamic limit 
N -> +00. In this limit P^ ~ 1/N and P^ - 1/N 2 . On 
the other hand, P x ~ 1 and |F(i j)\ ~ G ~ 1/JV (see 
Appendix IA1. 

The term in the l.h.s. of Eq. (fTO)) is of order 1, and 
the term in the r.h.s. is of order 1/N. Let us now con- 
sider the terms in Eq. (|11[) one by one. The first four 
terms correspond to the Liouville equation. The Liouvillian 
C = £0 + £12 + (£) describes the complete two-body prob- 
lem, including the inertial motion, the interaction between 
the particles (1, 2) and the mean field produced by the other 
particles. The terms £0 and (£) are of order 1/N while the 
term £ 12 is of order 1/N 2 . Therefore, the interaction term 
£i 2 can a priorF^I be neglected in the Liouvillian. This cor- 
responds to the weak coupling approximation where only 
the mean field term £0 + (£) is retained. The fifth term in 
Eq. (jllj) is a source term S expressible in terms of the one- 
body distribution; it is of order 1/N. If we consider only the 
mean field Liouvillian £0 + (£) and the source term S, as 
we shall do in this paper, we can obtain a kinetic equation 
for stellar systems that is the counterpart of the Landau 
equation in plasma physics. The sixth term is of order 1/N 
and it corresponds to collective effects (i.e. dressing of the 
particles by the polarization cloud). In plasma physics, this 
term leads to the Lenard-Balescu equation. It takes into ac- 
count dynamical screening and regularizes the divergence 
at large scales that appears in the Landau equation. In the 
case of stellar systems, there is no large-scale divergence be- 
cause of the spatial inhomogeneity of the system. Therefore, 
collective effects are less crucial in the kinetic theory of 
stellar systems than in plasma physics. However, this term 
has been properly taken into account by Heyvaerts (2010) 
who obtained a kinetic equation of stellar systems that is 
the counterpart of the Lenard-Balescu equation in plasma 
physics. The last two terms are of the order 1/A'' 2 and they 
will be neglected. In particular, the three-body correlation 
function P3, of order 1/N 2 , can be neglected at the order 
1/A. In this way, the hierarchy of equations is closed and 



Actually, the interaction term becomes large at small scales 
so its effect is not totally negligible (in other words, the expan- 
sion in terms of 1/JV is not uniformly convergent). In particular, 
the interaction term £12 must be taken into account in order to 
describe strong collisions with small impact parameter that lead 
to large deflections very different from the mean field trajectory 
corresponding to £0 + (£). We shall return to this problem in 
Sec. 1231 
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a kinetic equation involving only two-body encounters can 
be obtained. 

If we introduce the notations / = NmP\ (distribution 
function) and g = N 2 P' 2 (two-body correlation function), 
we get at the order 1/N: 

^■(1)+vi~(1) + — 
at or i 



J F(2->l)ff(l,2)dx a , 



(12) 



ldg 
2 dt 

^F(2 



(1,2) 



■vi -^(1,2)+ (F)(1)- ^(1,2) 



1 



Of 



F(3 



l) ff (2, 3) dxa 



(l)/(2) + F(2 
9/ 



1) 



9vi 



9vi 

(l) + (lo2) 



(1,2) 
0. 



(13) 



We have introduced the mean force (by unit of mass) cre- 
ated on star 1 by all the other stars 



(F)(1) = J F(2 



1)M, X2 



-V*(l), 



(14) 



and the fluctuating force (by unit of mass) created by star 
2 on star 1: 



F(2 



1) 



F(2 



1) 



iv (F)(1) - 



(15) 



Equations (|12[) and (fTB")) are exact at the order 1/N. They 
form the right basis to develop the kinetic theory of stellar 
systems at this order of approximation. Since the collision 
term in the r.h.s. of Eq. ([T2"j) is of order 1/N, we expect 
that the relaxation time of stellar systems scales as ~ Ntn 
where in is the dynamical time. As we shall see, the discus- 
sion is more complicated due to the presence of logarithmic 
corrections in the relaxation time and the absence of strict 
statistical equilibrium state. 



+oo: the Vlasov equation (collision/ess 



2.3. The limit N 
regime ) 

In the limit N — > +oo, for a fixed interval of time [0,T] 
(any), the correlations between stars can be neglected. 
Therefore, the mean field approximation becomes exact and 
the iV-body distribution function factorizes in N one-body 
distribution functions 



TV 



P N (xi,...,x N ,t) = JJPi(xj,i). 



(16) 



Substituting the factorization (|16[) in the Liouville equation 
(|2|), and integrating over X2, X3, xjy, we find that the 
smooth distribution function /(r, v, i) = NmPi(r,v,t) is 
the solution of the Vlasov equation 



dt 



2L 

Or 



-V$, A$ = AttG / fdv 



(17) 



This equation also results from Eq. (|12j) if we neglect the 
correlation function g(l, 2) in the r.h.s. and replace N — 1 
byN. 



The Vlasov equation describes the collisionless evolu- 
tion of stellar systems for times shorter than the relaxation 
time ~ Ntu- In practice N ^> 1 so that the domain of va- 
lidity of the Vlasov equation is huge (see the end of Sec. 
12. 8[) . As recalled in the Introduction, the Vlasov-Poisson 
system develops a process of phase mixing and violent re- 
laxation leading to a quasi-stationary state (QSS) on a very 
short timescale, of the order of a few dynamical times to- 
Elliptical galaxies are in such QSSs. Lyndcn-Bcll (1967) de- 
veloped a statistical mechanics of the Vlasov equation in or- 
der to describe this process of violent relaxation and predict 
the QSS achieved by the system. Unfortunately, the pre- 
dictions of his statistical theory are limited by the problem 
of incomplete relaxation. Kinetic theories of violent relax- 
ation, which may account for incomplete relaxation, have 
been developed by Kadomtscv and Pogutsc (1970), Scverne 
and Luwel (1980), and Chavanis (1998,2008). 

2.4. The order 0(1 /N): the generalized Landau equation 
(collisional regime) 

If we neglect strong collisions and collective effects, the first 
two equations (|12p and (fT3")) of the BBGKY hierarchy re- 
duce to 



f(0>+v.g<0> 



^V><°>!<0> 



d 



F(l->0) 5 (0,l)d Xl , (18) 



v>,i: 

or 
-V(l 



(F>(0). ^(0,1) + (0^1) 
0)-^(0)/(1) + (0 1). (19) 



The first equation gives the evolution of the one-body dis- 
tribution function. The l.h.s. corresponds to the (Vlasov) 
advection term. The r.h.s. takes into account correlations 
(finite V effects, graininess, discreteness effects) between 
stars that develop due to their interactions. These correla- 
tions correspond to "collisions" . 

Equation (|19|) may be viewed as a linear first order dif- 
ferential equation in time. It can be symbolically written 
as 



dg 
dt 



Cg = S[f], 



(20) 



where C = Cq + (C) is a mean field Liouvillian and S[f] 
is a source term S expressible in terms of the one-body 
distribution. This equation can be solved by the method of 
characteristics. Introducing the Green function 



G(t,t') =exp|- J £(r)dr| 



(21) 



constructed with the mean field Liouvillian C, we obtain 



ff (x,xi,f) = ~ / drG(t,t-T) 
™ Jo 

P(l 0). A + P(0 ->!)._!_ 

x/(x,f-r)/(xi,f-r), 



(22) 
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where we have assumed that no correlation is present ini- 
tially so that g(x, xi, £ = 0) = (if correlations are present 
initially, it can be shown that they are rapidly washed out). 
Substituting Eq. $FZ§ in Eq. dTHJ) , we obtain 



dt 



df N - 1 



Or 



N 



df 
dv 



■^JdrJ dr 1 dv 1 F"(1^0)G{t,t-r) 



F V {1 0) 



dv u 



F v (0 



dv\ 



/(r,v,£-r) — (ri,vi,t-r) 
m 



(23) 



In writing this equation, we have adopted a Lagrangian 
point of view. The coordinates appearing after the 
Greenian must be viewed as (£ — r) = (£) — J Q ds Vi (t — 

s) ds and v^(£ — r) = Vi(t) — / Q T ds (F)(rv(£ — s),t — s) ds. 
Therefore, in order to evaluate the integral (|2"3"]l . we must 
move the stars following the trajectories determined by the 
self-consistent mean field. 

The kinetic equation (f2"3"]l is valid at the order 1 /N so it 
describes the "collisional" evolution of the system (ignoring 
collective effects) on a timescale of order Nto- Equation 
(|2"3"|) is a non-Markovian integro-diffcrential equation. It 
takes into account derealizations in space and time (i.e. 
spatial inhomogencity and memory effects). Actually, the 
Markovian approximation is justified in the N — > +oo limit 
because the timescale ~ Nto over which /(r,v,£) changes 
is long compared to the correlation time r corr ~ over 
which the integrand in Eq. (|2"3")l has significant supporQ- 
Therefore, we can compute the correlation function (|2"2"T) by 
assuming that the distribution function is "frozen" at time 
£. This corresponds to the Bogoliubov ansatz in plasma 
physics. If we replace /(r, v,£ — r) and /(ri,Vi,f — r) by 
/(r, v,£) and /(ri,Vi,£) in Eq. (|23|) and extend the time 
integral to +oo, we obtain 

df df JV-1.„. df 

— + v • — H (F) • — 

dt dr N dv 

d f +ca f 

= W ^y^ v ^(^°) 



xG(£,£-r) 



F v (l -> 0)— + #"(0 -> I) — 
' dv v v y dv\ 



/(r,v,£) — (n,vi,i). 
m 



(24) 



Similarly, we can compute the trajectories of the stars by 
assuming that the mean field is independent on r and equal 
to its value at time £ so that r j (£ — r) = r, (£) — J Q (is Vj (£ — 
s) ds and v, (£ — t) = Vj (£ ) — / T (F) (r, (£ — s) , £) ds. 

The structure of the kinetic equation (f2~4")l has a clear 
physical meaning. The l.h.s. corresponds to the Vlasov ad- 
vection term due to mean field effects. The r.h.s. can be 



11 It is sometimes argued that the Markovian approxima- 
tion is not justified for stellar systems because the force 
auto-correlation decreases slowly as 1/t (Chandrasekhar 1944). 
However, this result is only true for an infinite homogeneous 
system (see Appendix [Gj . For spatially inhomogeneous distri- 
butions, the correlation function decreases more rapidly and 
the Markovian approximation is justified (Severne and Haggerty 
1976). 



viewed as a collision operator Cm [/] taking finite N effects 
into account. For N — > +oo, it vanishes and we recover the 
Vlasov equation. For finite N, it describes the cumulative 
effect of binary collisions between stars. The collision oper- 
ator is a sum of two terms: A diffusion term and a friction 
term. The coefficients of diffusion and friction are given by 
generalized Kubo formulae, i.e. they involve the time inte- 
gral of the auto-correlation function of the fluctuating force. 
The kinetic equation (|2"4"|) bears some resemblance with the 
Fokkcr-Planck equation. However, it is more complicated 
since it is an integro- differential equation, not a differential 
equation (see Section [3]). 

Equation ([M]) may be viewed as a generalized Landau 
equation. Since the spatial inhomogeneity and the finite 
extension of the system are properly taken into account, 
there is no divergence at large scales. There is, however, 
a logarithmic divergence at small scales which is due to 
the neglect of the interaction term £12 in the Liouvillian 
(see Sec. 12. 6[) . At large scales (i.e. for large impact param- 
eters), this term can be neglected and the trajectories of 
the stars are essentially due to the mean field. Thus, we 
can make the weak coupling approximation leading to the 
Landau equation. This approximation describes weak col- 
lisions. However, at small scales (i.e. for small impact pa- 
rameters), we cannot ignore the interaction term £12 in the 
Liouvillian anymore and we have to solve the classical two- 
body problem. This is necessary to correctly describe strong 
collisions for which the trajectory of the particles deviates 
strongly from the mean field motion. When the mean field 
Greenian G (constructed with £ + (£)) is replaced by the 
total Greenian G' (constructed with £ + £i2 + (£)), taking 
into account the interaction term, the generalized Landau 
equation (|24[) is replaced by a more complicated equation 
which can be viewed as a generalized Boltzmann equation. 
This equation is free of divergence (since it takes both spa- 
tial inhomogeneity and strong collisions into account) but 
it is unnecessarily complicated because it does not exploit 
the dominance of weak collisions over (rare) strong colli- 
sions for the gravitational potential. Indeed, a star suffers a 
large number of weak distant encounters and very few close 
encounters. A better practical procedure is to use the gen- 
eralized Landau equation (f2"4"]) with a cut-off at small scales 
in order to take into account our inability to describe strong 
collisions by this approach. 

The generalized Landau equation (|24[) was derived by 
Kandrup (1981) from the Liouvillc equation by using the 
projection operator formalism. It can also be derived from 
the Liouville equation by using the BBGKY hierarchy or 
from the Klimontovich equation by making a quasilinear 
approximation (Chavanis 2008). 

2.5. The Vlasov-Landau equation 

Self-gravitating systems are intrinsically spatially inhomo- 
geneous. However, the collision operator at position r can 
be simplified by making a local approximation and perform- 
ing the integrals as if the system were spatially homoge- 
neous with the density p = p(r). This amounts to replacing 
/(ri,vi,£) by /(r,Vi,£) and F{i -> 3) by F(i -> j) in Eq. 
(|24[). This local approximation is motivated by the work 
of Chandrasekhar and von Neumann (1942) who showed 
that the distribution of the gravitational force is a Levy 
law (called the Holtzmark distribution) dominated by the 
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contribution of the nearest neighbor. With this local ap- 
proximation, Eq. ([24)) becomes 

df df N-l^, df 

— + v • — H (F) ■ — 

dt dr N dv 

= Ho dT / dridviF/i ( 1 ^°) G °(*'*- r ) 
^(l-^0)(A_^) /(riVit) Z (rjVlit)i ( 25 ) 

where we have used F(0 -> 1) = — F(l -> 0). The Greenian 
Go corresponds to the free motion of the particles associ- 



ated with the Liouvillian Co. Using Eqs. (jC.ip and (IC.2I) 

the foregoing equation can be rewritten as 

df df N-l.„, df 

- r- V - ~ 1 (F) • — 

dt 3r JV u 3v 



O.t-r) 



Xl ^-^)^ V '*)^ r ' V1 '^ (26) 

where F(l — > 0, t — r) is expressed in terms of the 
Lagrangian coordinates. The integrals over r and ri can 
be calculated explicitly (see Appendix [C]). We then find 
that the evolution of the distribution function is governed 
by the Vlasov-Landau equation 

df df N-l^, df , n 

x^T / m ( k -^k)( fl §L-f^ M , (27) 

where we have noted w = v — vi, / = f(r,v,t), f\ = 
f(r,Vi,t), and where (27r) 3 -u(fc) = —4irG/k 2 represents the 
Fourier transform of the gravitational potential. Under this 
form, we see that the collisional evolution of a stellar system 
is due to a condition of resonance kv = kv' (with v/V) 
encapsulated in the ^-function. This <5-function expresses 
the conservation of energy. 

The Vlasov-Landau equation can also be written as (see 
Appendix [C]) : 



dt dr 



K" v = A- 
where 



N dv 
d 



df t dh 



dv v dv\ 



2irmG 2 In A, 



ur 



(28) 
(29) 



interaction only appears in the constant A which merely de- 
termines the relaxation time. The structure of the Landau 
equation is independent on the potential. The Landau equa- 
tion was originally derived from the Boltzmann equation in 
the limit of weak deflections |Av| <C 1 (Landau 1936 j3- In 
the case of plasmas, the system is spatially homogeneous 
and the advection term is absent in Eq. (|2"5]l . In the case of 
stellar systems, when we make the local approximation, the 
spatial inhomogencity of the system is only retained in the 
advection term of Eq. (|2"5]l . This is why this kinetic equa- 
tion is referred to as the Vlasov-Landau equation. This is 
the fundamental kinetic equation of stellar systems. 



2.6. Heuristic regularization of the divergence 

To obtain the Vlasov-Landau equation ([28)1 , we have made 
a local approximation. This amounts to calculating the col- 
lision operator at each point as if the system were spatially 
homogeneous. As a result of this homogeneity assump- 
tion, a logarithmic divergence appears at large scales in the 
Coulombian factor (|30[) . In plasma physics, this divergence 
is cured by the Debye shielding. A charge is surrounded by 
a polarization cloud of opposite charges which reduces the 
range of the interaction. When collective effects are prop- 
erly taken into account, as in the Lenard-Balescu equation, 
no divergence appears as large scales and the Debye length 
arises naturally. Heuristically, we can use the Landau equa- 
tion and introduce an upper cut-off at the Debye length 
A_d = (ksT/ne 2 ) 1 ^ 2 . For self-gravitating systems, there is 
no shielding and the divergence is cured by the finite extent 
of the system. The interaction between two stars is only lim- 
ited by the size of the system. When spatial inhomogencity 
is taken into account, as in the generalized Landau equa- 
tion ([2"4")h no divergence occurs at large scales. Heuristically, 
we can use the Vlasov-Landau equation ([28)) and introduce 
an upper cut-off at the Jeans length Aj ~ (/csT/Gm 2 ™) 1 / 2 
which gives an estimate of the system size R. 

The Coulombian factor ([3U)) also diverges at small 
scales. As explained previously, this is due to the neglect 
of strong collisions that produce important deflections. 
Indeed, for collisions with low impact parameter, the mean 
field approximation is clearly irrelevant and it is necessary 
to solve the two-body problem exactly. Heuristically, we 
can use the Landau equation and introduce a lower cut- 
off at the gravitational Landau length Xl ~ Gm/v^ ~ 
Gm 2 /(fcsT) (the gravitational analogue of the Landau 
length Al ~ e 2 /mv 2 a ~ e 2 /fc^T in plasma physics) which 
corresponds to the impact parameter leading to a deflection 
at 90°. 

Introducing a large-scale cut-off at the Jeans length Aj 
and a small-scale cut-off at the Landau length A l , and not- 
ing that Al ~ l/(nAj), we find that the Coulombian fac- 



lnA = 



dk/k, 



(30) 



is the Coulomb factor that has to be regularized with 
appropriate cut-offs (see Section 12. 6[) . The r.h.s. of Eq. 
(J25J is the original form of the collision operator given by 
Landau (1936) for the Coulombian interactiorP^I. It applies 
to weakly coupled plasmas. We note that the potential of 



12 The case of plasmas is recovered by making the substitution 
Gm 2 o -e 2 leading to K = (2ire 4 /m 3 ) In A. 



13 We note that Eq. ((25)) with Go replaced by the total 
Greenian G' taking into account the interaction term is equiv- 
alent to the Boltzmann equation. Indeed, for spatially homo- 
geneous systems, the Boltzmann equation can be derived from 
the BBGKY hierarchy (|10|) and Hll|) by keeping the Liouvillian 
C — Co + C12 describing the two-body problem exactly and 
the source S (Balescu 2000). Therefore, the procedure used by 
Landau which amounts to expanding the Boltzmann equation 
in the limit of weak deflections is equivalent to the one presented 
here that starts from the BBGKY hierarchy and neglects £12 in 
the Liouvillian. 
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tor scales as InA ~ ln(Aj/Ai) ~ ln(nAj) ~ In AT where 



N ■ 



R 3 is the number of stars in the clustciP^I. 



2. 7. Properties of the Vlasov- Landau equation 

The Vlasov-Landau equation conserves the mass M = 

J f drdv and the energy E = J f^drdv + \\ p&dr. 
It also monotonically increases the Boltzmann entropy 
S = - J(f /m)ln(f /m) drdv in the sense that S > (H- 
theorem). Due to the local approximation, the proof of 
these properties is the samcFl as for the spatially homo- 
geneous Landau equation (Balescu 2000). Because of these 
properties, we might expect that a stellar system will relax 
towards the Boltzmann distribution which maximizes the 
entropy at fixed mass and energy. However, we know that 
there is no maximum entropy state for an unbounded self- 
gravitating system (the Boltzmann distribution has infinite 
mass). Therefore, the Vlasov-Landau equation does not re- 
lax towards a steady state and the entropy does not reach 
a stationary value. Actually, the entropy increases perma- 
nently as the system evaporates. But since evaporation is 
a slow process, the system may achieve a quasistationary 
state that is close to the Boltzmann distribution. A typical 
quasistationary distribution is the Michic-King model 



(31) 



where e = v 2 /2 + <!>(r) is the energy and j=rxv the an- 
gular momentum. This distribution takes into account the 
escape of high energy stars and the anisotropy of the ve- 
locity distribution. It can be derived, under some approx- 
imations, from the Vlasov-Landau equation by using the 
condition that / = if the energy of the star is larger 
than the escape energy e m (Michie 1963, King 1965). The 
Michie-King distribution reduces to the isothermal distri- 
bution / cx e^' 3me for low energies. In this sense, we can 
define a "relaxation" time for a stellar system. From the 
Vlasov-Landau equation (f2"8")) , we find that the relaxation 
time scales as 



t R 



i 2 G 2 In AT' 



(32) 



where v m is the root mean square (r.m.s.) velocity. 
Introducing the dynamical time tjj ~ Xj/v m ~ R/v m , we 
obtain the scaling 



N 



(33) 



The fact that the ratio between the relaxation time and the 
dynamical time depends only on the number of stars and 
scales as N/ In N was first noted by Chandrasckhar (1942). 

A simple estimate of the evaporation time gives t evap ~ 
136ifl (Ambartsumian 1938, Spitzer 1940). More precise 
values have been obtained by studying the evaporation 
process in an artificially uniform cluster (Chandrasekhar 
1943b, Spitzer and Harm 1958, Michie 1963, King 1965, 



Similarly, in plasma physics, introducing a large-scale cut-off 
at the Debye length \d and a small-scale cut-off at the Landau 
length Al, and noting that Al ~ l/(n\ 2 D ), we find that the 
Coulombian factor scales as InA ~ ln(A_o / 'Al) ~ ln(nA|j) where 
A ~ n\% is the number of electrons in the Debye sphere. 
15 We recall that the Vlasov advection term conserves mass, 
energy, and entropy. 



Lemou and Chavanis 2010) or in a more realistic inhomo- 
geneous cluster (Henon 1961). Since t evap 3> tn, we can 
consider that the system relaxes towards a steady distribu- 
tion of the form (|3"Tj) on a timcscale tu and that this distri- 
bution slowly evolves on a longer timescale as the stars es- 
capj^l. The characteristic time in which the system's stars 
evaporate is t evap . The evaporation is one reason for the 
evolution of stellar systems. However, as demonstrated by 
Antonov (1962) and Lynden-Bell and Wood (1968), stel- 
lar systems may evolve more rapidly as a result of the 
gravothermal catastrophe. In that case, the Michie-King 
distribution changes significatively due to core collapse. 
This evolution has been described by Cohn (1980), and 
it leads ultimately to the formation of a binary star sur- 
rounded by a hot halo. Such a configuration can have ar- 
bitrarily large entropy. Cohn (1980) finds that the entropy 
increases permanently during core collapse, confirming that 
the Vlasov-Landau equation has no equilibrium state. 

Even if the system were confined within a small box so 
as to prevent both the evaporation and the gravothermal 
catastrophe, there would be no statistical equilibrium state 
in a strict sense because there is no global entropy maxi- 
mum (Antonov 1962). A configuration in which some sub- 
set of the particles are tightly bound together (e.g. a binary 
star), and in which the rest of the particles shares the en- 
ergy thereby released, may have an arbitrary large entropy. 
However, such configurations, which require strong corre- 
lations, are generally reached very slowly (on a timescale 
much larger than (AT/lnAT)^) due to encounters involv- 
ing many particles. To describe these configurations, one 
would have to take high order correlations into account in 
the kinetic theory. These configurations may be relevant 
in systems with a small number of stars (Chabanol et al. 
2000) but when N is large the picture is different. On a 
timcscale of the order of (N/ In N)tr> the one-body distri- 
bution function is expected to reach the Boltzmann dis- 
tribution which is a local entropy maximum. This state is 
"metastable" but its lifetime is expected to be very large, 
scaling as e N , so that it is stable in practice (Chavanis 
2006). In this sense, there exist "true" statistical equilib- 
rium states for self-gravitating systems confined within a 
small box. However, we may argue that this situation is 
highly artificial. 

2.8. Dynamical evolution of stellar systems: a short review 

Using the kinetic theory, we can identify different phases in 
the dynamical evolution of stellar systems. 

A self-gravitating system initially out-of-mechanical 
equilibrium undergoes a process of violent collisionless re- 
laxation towards a virialized state. In this regime, the dy- 



10 The distribution (|31|) is not steady in the sense that the co- 
efficients A, p, and r a slowly vary in time as the stars escape 
and the system loses mass and energy. However, the distribu- 
tion keeps the same form during the evaporation process. Henon 
(1961) found a self-similar solution of the orbit-averaged- Fokker- 
Planck equation. He showed that the core contracts as the clus- 
ter evaporates and that the central density is infinite (the struc- 
ture of the core resembles the singular isothermal sphere with a 
density p oc r~ 2 ). He argued that the concentration of energy, 
without concentration of mass, at the center of the system is 
due to the formation of tight binary stars. The invariant profile 
found by Henon does not exactly coincide with the Michie-King 
distribution (introduced later), but is reasonably close. 
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namical evolution of the cluster is described by the Vlasov- 
Poisson system. The phenomenology of violent relaxation 
has been described by Henon (1964), King (1966), and 
Lyndcn-Bcll (1967). Numerical simulations that start from 
cold and clumpy initial conditions generate a quasi sta- 
tionary state (QSS) that fits the de Vaucouleurs i? 1 / 4 law 
for the surface brightness of elliptical galaxies quite well 
(van Albada 1982). The inner core is almost isothermal 
(as predicted by Lynden-Bell 1967) while the velocity dis- 
tribution in the envelope is radially anisotropic and the 
density profile decreases as r~ 4 . One success of Lynden- 
Bell's statistical theory of violent relaxation is to explain 
the isothermal core of elliptical galaxies without recourse 
to "collisions" . By contrast, the structure of the halo can- 
not be explained by Lynden-Bell's theory as it results from 
an incomplete relaxation. Models of incompletely relaxed 
stellar systems have been elaborated by Bertin and Stiavelli 
(1984), Stiavelli and Bertin (1987), and Hjorth and Madsen 
(1991). These theoretical models nicely reproduce the re- 
sults of observations or numerical simulations (Londrillo et 
al. 1991, Trenti et al. 2005). In the simulations, the initial 
condition needs to be sufficiently clumpy and cold to gen- 
erate enough mixing required for a successful application of 
the statistical theory of violent relaxation. Numerical sim- 
ulations starting from homogeneous spheres (see, e.g., Roy 
and Perez 2004, Levin et al. 2008, Joyce et al. 2009) show 
little angular momentum mixing and lead to different re- 
sults. In particular, they display a larger amount of mass 
loss (evaporation) than simulations starting from clumpy 
initial conditions. Clumps thus help the system to reach a 
"universal" final state from a variety of initial conditions, 
which can explain the similarity of the density profiles ob- 
served in elliptical galaxies. 

On longer timescales, encounters between stars must be 
taken into account and the dynamical evolution of the clus- 
ter is governed by the Vlasov-Landau-Poisson system. The 
first stage of the collisional evolution is driven by evap- 
oration. Due to a series of weak encounters, the energy 
of a star can gradually increase until it reaches the local 
escape energy; in that case, the star leaves the system. 
Numerical simulations (Spitzer 1987, Binney and Tremaine 
2008) show that during this regime the system reaches a 
quasi-stationary state that slowly evolves in amplitude due 
to evaporation as the system loses mass and energy. This 
quasi stationary distribution function is close to the Michie- 
King model (|3lTl . The system has a core- halo structure. The 
core is isothermal while the stars in the outer halo move 
in predominantly radial orbits. Therefore, the distribution 
function in the halo is anisotropic. The density follows the 
isothermal law p ~ r~ 2 in the central region (with a core of 
almost uniform density) and decreases as p ~ r _7/2 in the 
halo (for an isolated cluster with e m = 0). Due to evapo- 
ration, the halo expands while the core shrinks as required 
by energy conservation. During the evaporation process, the 
central density increases permanently. At some point of the 
evolution, the system unde rgo es an instability related to the 
Antonov (1962) instabilhVsrl and the gravothermal catas- 
trophe sets in (Lynden-Bell and Wood 1968). This instabil- 



17 The collisional evolution of the system can be measured pre- 
cisely in terms of the scale escape energy xo = 3$ (0) (0) . The 
onset of instability corresponds to xo — 9.3 which is the value 
at which the King model becomes thermodynamically unstable 
(Cohn 1980). 



ity is due to the negative specific heat of the inner system 
that evolves by losing energy and thereby growing hotter. 
The energy lost is transferred outward by stellar encoun- 
ters. Hence the temperature always decreases outward, and 
the center continually loses energy, shrinks, and heats up. 
This leads to core collapse. Mathematically speaking, core 
collapse would generate a finite time singularity. When the 
evolution is modeled by the orbit-averagcd-Fokkcr-Planck 
equation, Cohn (1980) finds that the collapse is self-similar, 
that the central density becomes infinite in a finite time, 
and that the density behaves as p ~ r - 2 -23^ invariant 
profile found by Cohn differs from the Michie-King distribu- 
tion (for which p ~ r~ 2 ) beyond a radius of about 10r core . 
Larson (1970) and Lyndcn-Bcll and Egglcton (1980) find 
similar results by modeling the evolution of the system by 
fluid equations. Alternatively, when the evolution is mod- 
eled by the original Vlasov-Landau equation, Lancellotti 
and Kiessling (2001) argue that the density should behave 
a,s p oc r~ 3 in the final stage of the collapse. In all cases, the 
authors find a singular density profile at the collapse time 
that is integrable (or diverges logarithmically) at r = 0. 
This means that the core contains very little mass. In re- 
ality, if we come back to the iV-body system, there is no 
singularity and core collapse is arrested by the formation 
of binary stars due to three-body collisions. These binaries 
can release sufficient energy to stop the collapse and even 
drive a re-expansion of the cluster in a post-collapse regime 
(Inagaki and Lynden-Bell 1983). Then, in principle, a series 
of gravothermal oscillations should follow (Bettwieser and 
Sugimoto 1984). 

At the present epoch, small groups of stars such as 
globular clusters (N ~ 10 5 , t D — 10 5 yr, age ~ 10 10 yr, 
tn ~ 10 10 yr) are in the collisional regime. They are ei- 
ther in quasistationary states described by the Michie-King 
model or experiencing core collapse. By contrast, large clus- 
ters of stars like elliptical galaxies (N ~ 10 11 , tn ~ 10 s yr, 
age ~ 10 10 yr, tn ~ 10 19 yr) are still in the collisionless 
regime and their apparent organization is a result of an 
incomplete violent relaxation. 

3. Test star in a thermal bath 

3.1. The Fokker-Planck equation 

We now consider the relaxation of a "test" star (tagged 
particle) evolving in a steady distribution of "field" stariQ. 
Due to the encounters with the field stars, the test star has a 
stochastic motion. We call P(r, v, t) the probability density 
of finding the test star at position r with velocity v at 
time t. The evolution of P(r, v, t) can be obtained from the 
generalized Landau equation ([24]) by considering that the 



In plasma physics, this steady distribution is the Boltzmann 
distribution of statistical equilibrium which is the steady state 
of the Landau equation. In the case of stellar systems, there is 
an intrinsic difficulty since no statistical equilibrium state exists 
in a strict sense: The Boltzmann distribution has infinite mass, 
and the Vlasov-Landau equation has no steady state because of 
the escape of high energy stars. However, we have seen that the 
system can reach a quasi steady distribution (e.g. a Michie-King 
distribution) and that this distribution changes on an evapora- 
tion timescale that is long with respect to the collisional relax- 
ation time. Therefore, we can consider that this distribution is 
steady on the collisional timescale over which the Fokker-Planck 
approach applies. 
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distribution function of the field stars is fixed. Therefore, we 
replace /(r,v,i) by P(r,v,i) and /(n,vi,i) by /(ri,vi) 
where /(ri , vi) is the steady distribution of the field starj^l. 
This procedure transforms the intcgro-differcntial equation 
((Ml) into the differential equation 

dP dP N - 1 DP 
~dt +V '!fr + N '"for 



d 



+00 



xG(t,t-r) 



(It 
8 



N ' 
dridviP"(l 



0) 



^(1^0)— +F^(0^1) 



f 



P(r,v,i)-(ri, Vl ) 



where (F)(r) = — V$(r) is the static mean force created by 
the field stars with density p(r\). This equation does not 
present any divergence at large scales. 

If wc make a local approximation and use the Vlasov- 
Landau equation (|26[) . wc obtain 



dP 
~dt 



dP 
dr 



N - 1 



_d_ 



/■+00 /■ 
J dr J drxdviF^il 



_d_ 



d 



N 

0,t)F v (l - 
/ 



dP 
c/v 

0,t-r) 



P(r,v,t)^(r, Vl ), (35) 



Denoting the advection operator by d/dt, Eq. (|35|) can be 
written in the form of a Fokkcr-Planck equation 



dt dv^ \ dv u 
involving a diffusion tensor 



rr pol 



(36) 



DtiV = _ 
m 



+00 



dr 



xP(l 

and a friction force 

f+00 



d ri dviP^(l -> 0,t) 
->0,t-T)/(r,vi), 



(37) 



p; oi = l^ 

xK(l -> 0,t 



0,t) 



r) cK (r ' Vl) ' 



(38) 



If we directly start from Eq. ([27)1 . which amounts to 
performing the integrals over r and ri in the previous ex- 
pressions, we obtain the Fokker-Planck equation 

dP dP N - 1 

~dt +v 



d 
dvf* 



dr 

k»k v 5(k • w)ti(fc) 2 



N 
fi 



— = 7r(27rrm 



<9p 



p 



dfi 



dvjdk, (39) 



We can understand this procedure as follows. Equations (|24l) 
and (|34p govern the evolution of the distribution function of 
a test star (described by the coordinates r and v) interact- 
ing with field stars (described by the running coordinates ri 
and vi). In Eq. 1)2411 . all the stars are equivalent so the dis- 
tribution of the field stars /(ri,vi,<) changes with time ex- 
actly like the distribution of the test star /(r, v,£). In Eq. 1)340. 
the test star and the field stars are not equivalent since the 
field stars form a "bath". The field stars have a steady (given) 
distribution /(ri,vi) while the distribution of the test star 
/(r, v,t) = NmP(r,v,t) changes with time. 



with the diffusion and friction coefficients 

= 7r(27r) 3 m J k»k v 8{k • w)u(fc) 2 /i dv x dk, (40) 



F£, = 7 r(2 7 r) 3 ?7i J k»k v 6{k ■ w)u(k) 2 ^ d^dk. (41) 

The diffusion tensor D^ v results from the fluctuations 
of the gravitational force due to the granularities in the 
distribution of the field stars. It can be derived directly 
from the formula (see Appendix [G| : 



(34) D"" 



+ DC 



(F»{t)F v {t-T))dT, 



(42) 



deduced from Eq. ([44l -a). The friction force F po ; results 
from the retroaction of the field stars to the perturbation 
caused by the test star like in a polarization process. It 
can be derived from a linear response theory (Marochnik 
1968, Kalnajs 1971b, Kandrup 1983, Chavanis 2008). It 
will be called the "friction by polarization" to distinguish 
it from the total friction (see below). Eqs. (j35)) - ([4lT) have 
been derived within the local approximation. More general 
formulae, valid for fully inhomogeneous stellar systems, are 
given in Kandrup (1983) and Chavanis (2008). The friction 
force has also been calculated by Tremaine & Weinberg 
(1984), Bekenstein & Maoz (1992), Maoz (1993), Nelson & 
Tremaine (1999) using different approaches. 

Since the diffusion tensor depends on the velocity v of 
the test star, it is useful to rewrite Eq. (jSHJ) in a form that 
is fully consistent with the general Fokker-Planck equation 



dP _ 

dt dv^dv 
with 

DlLV = (Av»Av») 



d 2 v d „ 

-(D^ V P) - g-j:(PF frlcUon ), 



friction 



(Aw' 1 

2At ' * faction At 

By identification, we find that 
fa* = p» + 

friction pol 1 Qy 11 ' 



(43) 



(44) 



(45) 



Therefore, when the diffusion coefficient depends on the 
velocity, the total friction is different from the friction by 
polarization. Substituting Eqs. (00]) and flUJ) in Eq. (|35)l. 
and using an integration by parts, we find that the diffusion 
and friction coefficients can be written as 



(Av^Av"} 



2At 



7r(27r) 3 m J k»k v 5{k- w)M(fc) 2 /irfv 1( ik, 



(Av 1 * 
At 



7t(27t) 3 to / k^k v f 



d d 
xS(k ■ w)u(k) 2 dvxdk 



(46) 



(47) 



These expressions can be obtained directly from the equa- 
tions of motion by expanding the trajectories of the stars 
in powers of 1/N in the limit N — > +oo (Chavanis 2008). 
We recall that Eqs. ((46)) and ((47)) display a logarithmic di- 
vergence a small and large scales that must be regularized 
by introducing proper cut-offs as explained in Section 12.61 
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In astrophysics, the diffusion and friction coefficients of a 
star were first calculated by Chandrasekhar (1943a) from 
a two-body encounters theory (see also Cohen et al. 1950, 
Gasiorowicz et al. 1956, and Roscnbluth et al. 1957). The 
expressions obtained by these authors arc different from 
those given above, but they are equivalent (see Section 
13. 5[) . In plasma physics, the diffusion and friction coeffi- 
cients of a charge were first calculated by Hubbard (1961a) 
who took collective effects into account, thereby eliminat- 
ing the divergence at large scales. When collective effects 
are neglected, his expressions reduce to Eqs. (|4*B")) and (|4"T)) . 
On the other hand, strong collisions have been taken into 
account by Chandrasekhar (1943a) in astrophysics and by 
Hubbard (1961b) in plasma physics. In that case, there is 
no divergence at small impact parameters in the diffusion 
and friction coefficients, and the Landau length appears 
naturally. 

The two forms (f3"B"]) and (j4"3"|) of the Fokker-Planck equa- 
tion have their own interest. The expression (|43[) where 
the diffusion coefficient is placed after the two derivatives 
d 2 (DP) involves the total friction force F friction and the 
expression (|36D where the diffusion coefficient is placed be- 
tween the derivatives dDdP isolates the part of the friction 
Fp i due to the polarization. Astrophysicists are used to the 
form ([4"3"]l . However, it is the form ([3"B"f that stems from the 
Landau equation (j2"6")l . We shall come back to this observa- 
tion in Section [3.51 

From Eqs. (|40|) and (|4"T1) . we easily obtain 



— r pol- 



dD u,y 

dv v 

Combining Eq. (gSJ) with Eq. (gSJ) , we get 



friction 



2F 



pol • 



(48) 



(49) 



Therefore, the friction force F f riction is equal to twice the 
friction by polarization F po i (for a test star with mass m in- 
teracting with field stars with mass m f , this factor two is re- 
placed by (m + m/)/m; see Appendix IP")). This explains the 
difference of factor 2 in the calculations of Chandrasekhar 
(1943a) who determined F friction and in the calculations 
of Kalnajs (1971b) and Kandrup (1983) who determined 



3.2. The Einstein relation 

In the central region of the system, the distribution of the 
field stars is close to the Maxwell-Boltzmann distribution 



/(n,vi) 



(50) 



where /3 = 1/T is the inverse temperature and p(ji) cx 
e -^m$(ri) ^ e density given by the Boltzmann law. 
Therefore, the field stars form a thermal bath. Substituting 
the identity 



= -/^TO/iVi, 

OVi 



(51) 



in Eq. (|4ip . using the S- function to replace k • vi by k • v, 
and comparing the resulting expression with Eq. (|40j) . we 
find that 



(52) 



The friction coefficient is given by an Einstein relation ex- 
pressing the fluctuation-dissipation thcorcrrPI We empha- 
size that the Einstein relation is valid for the friction force 
by polarization F po ;, not for the total friction F/riction 
(we do not have this subtlety in standard Brownian the- 
ory where the diffusion coefficient is constant). Using Eq. 
(|4"2"|) . we can rewrite Eq. (f5"2"j) in the form 



r pol 



-firm? 



dT(F"(t)F V (t-T)). 



(53) 



This relation is usually called the Kubo formula. More gen- 
eral expressions of the Kubo formula valid for fully inho- 
mogeneous stellar systems are given in Kandrup (1983) and 
Chavanis (2008). Using the Einstein relation, the Fokker- 
Planck equation (f3"6"|) takes the form 



dP 
dt 



d 



( dP 

D^(v) —+pmPv u 
\ov v 



(54) 



where the diffusion coefficient is given by Eq. (|4TJ)) with Eq. 
(|50p . This equation is similar to the Kramers equation in 
Brownian theory (Kramers 1940) except that the diffusion 
coefficient is a tensor and that it depends on the velocity 
of the test star. For an isotropic distribution function (e.g., 
the Maxwellian) , it can be put in the form 

1 v^v" 1 
D^ = (D \\--D ± ) — + -D ± 6" v , (55) 

where D\\ and D± are the diffusion coefficients in the direc- 
tions parallel and perpendicular to the velocity of the test 
star. The friction by polarization can then be written as 

F pol = -D\\Pmv. (56) 

The friction is proportional and opposite to the velocity 
of the test star, and the friction coefficient is given by the 
Einstein relation £ = Dn(3m. The total friction is 



friction 



= -2D\\f3mv 



(57) 



This is the Chandrasekhar dynamical friction. 

The steady state of the Fokker-Planck equation ([54]) is 
the Maxwell-Boltzmann distribution (|50p . Since the Fokker- 
Planck equation admits an ii-theorem for the Boltzmann 
free energy (Riskcn 1989), one can prove that P(r, v, t) con- 
verges towards the Maxwell-Boltzmann distribution ([501) 
for t — > +oo. In other words, the test star acquires the 
distribution of the field stars (thermalization). We recall, 
however, that for self-gravitating systems, these results 
cannot hold everywhere in the cluster since the Maxwell- 
Boltzmann distribution does not exist globally. 

If we assume that the distribution P(y,t) of the test 
star is isotropic, the Fokker-Planck equation becomes 



dP 
~dt 



d_ 
dv 



/3mPv 



(58) 



20 The collisional evolution of a star, under the effect of two- 
body encounters, can be understood as an interplay between 
two competing effects. The fluctuations of the gravitational field 
induce a diffusion in velocity space which tends to increase the 
speed of the star. This effect is counterbalanced by a friction 
(dissipation) which results in a systematic deceleration along the 
direction of motion. As emphasized by Chandrasekhar (1943a, 
1949) in his Brownian theory of stellar motion, the Einstein 
relation (|52p guarantees that the Maxwell distribution (|50l) is a 
steady state of the Fokker-Planck equation 1)3611 . 
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It can be obtained from an effective Langevin equation 



2£>,.(i;)R(f), (59) 



where R(i) is a Gaussian white noise satisfying (R(t)) = 
and (R^R^it')) = S^S(t - t'). Since the diffusion co- 
efficient depends on the velocity, the noise is multiplica- 
tive. The force acting on the star consists in two parts: 
a part derivable from a smoothed distribution of matter 
— V$ and a residual random part due to the fluctuations 
of that distribution. In turn, the random part can be de- 
scribed by a Gaussian white noise multiplied by a velocity- 
dependent factor plus a friction force equal to the friction 
by polarization F po ; = — Du(v)j3mv and a spurious drift 
F spurious = (l/2)(3D||/c?v) due to the multiplicative noise. 
If we neglect the velocity dependence of the diffusion coeffi- 
cient, we recover the results of Chandrasckhar (1943a, 1949) 
obtained in his Brownian theory. 

3.3. The dif Fusion tensor for isothermal systems 

When the velocity distribution of the field stars is given by 
the Maxwellian distribution ((50)) . the diffusion tensor (|40|) 
can be calculated as follows. If we introduce the represen- 
tation 



S(x) 



Atx uu 



2tt' 



(60) 



for the 6- function in Eq. (|4U|) . the diffusion tensor can be 
rewritten as 



£>"" = -(2tt) 6 
2 y ' 



dt / dk^r?}(fc)V kv 7(kt), (61) 



where / is the three-dimensional Fourier transform of the 
velocity distribution. This equation can be directly ob- 
tained from the formula (|4"2j) or (|3"T|) . This shows that 
the auxiliary integration variable t in Eq. (|6ip represents 
the time. For the Maxwellian distribution ([5U)) . /(kt) is a 
Gaussian. If we perform the integration over t (which is the 
one-dimensional Fourier transform of a Gaussian), we find 
that the diffusion tensor can be expressed as 



5 /{3m\ 1/2 f k^k v , 
D^ v = tt(2ttY [ — — \ pm j —j^— u ( 



\{k) 2 e- f)m{J ^- dk. 



(62) 



Alternatively, this expression can be obtained from Eq. (j4"0"|) 
by introducing a cartesian system of coordinates for vi with 
the z-axis taken along the direction of k, and performing 
the integration. With the notation k = k/fc, Eq. (|62p can 
be rewritten as 

= t:(2t:) 3 (^j ^ pm £°° k 3 u(k) 2 dk 



where 

G^(x) = [ k"k v e-^ 2 dk. 



(63) 



(64) 



We note that the potential of interaction only appears in a 
multiplicative constant that fixes the relaxation time (see 
below). Using Eq. (|C.7|) . we get 



D> 



3 y /2 2pmG 2 \nA 
2m 



QP .v 




(65) 



where In A is the Coulombian factor (|3U)) and v 2 n = 3/ (/3m) 
is the mean square velocity of the field stars. The diffusion 
tensor may be written as 



3t R \ V 2 v m , 



(66) 



where tR is the local relaxation time defined below in Eq. 
(|74l) . This relation emphasizes the scaling D ~ v^/tR. 

Introducing a spherical system of coordinates with the 
z-axis in the direction of x, we can write the normalized 
diffusion tensor in the form 



i T M T f i 



(67) 



where 

2tt 3 / 2 
G,| = G(x), 



with 



2tt 3 / 2 

G x = [erf(ar) - G(x)}, 



t 2 e~ l dt 



1 

2^ 



erf(x) — < -j=e x 

\/7T 



The error function is defined by 

2 f x 2 
erf (2) = I e _t dt. 
V 7 *" Jo 

We have the asymptotic behaviors 

f, 



Gi,(0) 



G\\(x) 



r 3/2 



2tt 3 / 2 



(69) 
(70) 

(71) 
(72) 



i J x 
We note that G^(x) ~ G\\(Q)8^ V when |x| -> 0. 

3.4. The relaxation time 

We can use the preceding results to estimate the relax- 
ation time of the velocity distribution of the test particle 
towards the Maxwellian distribution (thcrmalization) . If we 
set x = a/ /3m/2v, the Fokker-Planck equation ([54]) can be 
rewritten as 



dP _ 1 d 
dt ~ t R dx^ 



G^{x) 



dP 

dx» 



2Px v 



where tR is the local relaxation time 



tR 



1/2 



pmG 2 In A 



(73) 



(74) 



The prefactor is equal to 0.482 (of course this numerical 
factor may vary depending on the definition of the relax- 
ation time). The relaxation time is inversely proportional 
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Fig. 1. Normalized diffusion coefficients G\\(x), G±(x) and 
friction force xG|| (x) for a thermal bath. The friction is 
maximum for x ~ 0.97, i.e. when the velocity of the test 
star is approximately equal to the r.m.s. velocity of the field 
stars. 



to the local density p{r). Therefore, the relaxation time is 
smaller in regions of high density (core) and larger in re- 
gions of low density (halo) . Introducing the dynamical time 
to = Xj/v m , we get 



tR 



N 



t D . 



(75) 



We note that the relaxation time of a test particle in a bath 
is of the same order as the relaxation time of the system as 
a whole (see Sec. I2.7[) . This property is not true anymore 
in one dimension (Eldridgc and Fcix 1962). 

We can also get an estimate of the relaxation time by 
the following argument (Spitzcr 1987). If the diffusion coef- 
ficient were constant, the typical velocity of the test star (in 
one spatial direction) would increase as ((Av) 2 )/3 ~ 2D»t. 
The relaxation time t r is the typical time at which the 
typical velocity of the test star has reached its equilibrium 
value (v 2 ){+oo) = 3/(m/3) = v 2 n so that ((Av) 2 )(t r ) = 
(v 2 )(+oo). Since D\\ depends on v, the description of the 
diffusion is more complex. However, the formula t r = 
v^/l&D^ivrn)] resulting from the previous arguments with 
Z?|| = D\\{v m ) should provide a good estimate of the re- 
laxation time. Using Eq. (|65p and comparing with Eq. 
1(74) we obtain t r = K 3 t R , where K 3 = \/[AG\\{^Jzj2)}. 
Numerically, K 3 = 0.13587547.... 

Finally, we can estimate the relaxation time by t' r = 
where £ is the friction coefficient. Using the Einstein 
relation £ = 2D\\f3m [sec Eq. ((57)) ] with D\\ = D\\(v m ) we 
find that t' = t r . 



3.5. The Rosenbluth potentials 

It is possible to obtain simple expressions of the diffusion 
and friction coefficients for any isotropic distribution of the 
bath. If wc start from the expression (|28|) of the Vlasov- 
Landau equation, we find that the Fokker-Planck equation 
(f3"5() can be written as 



dP 
~dt 



d 



n dv v dv\ 



dv\. 



(76) 



igeneous stellar systems without collective effects 
The diffusion and friction coefficients are given by 



= / K^h dvi=A h 



dv u (77) 



friction 



J dv 
Using the identities 

2 

-((' 



f w^ L 

/idvi = -AA J /ydvi. (78) 



= A 

and 
dK»" 



d 2„ 



dv^dv v 



(79) 



-2A— = 2A— - 



(80) 

ov>* \w ) 

the coefficients of diffusion and friction can be rewritten 
d 2 9 



as 



dv^dv 



dh 

F friction = 2F po ; = 4A—(y), 

where 



(81) 
(82) 



5(v) = J /(vi)|v-vi|dvi, h(v) = J J^-dv!, (83) 

are the so-called Rosenbluth potentials (Rosenbluth et al. 
1957). 

If the field particles have an isotropic velocity distribu- 
tion, the Rosenbluth potentials take the particularly simple 
form (see, e.g., Binney and Tremaine 2008): 



h(y) = 47T 



2 rv r + co 

- / f{v\)vldv\ + / f(vi)vidvi 



(84) 



K«) = it 



J ( — + vvij f(vi)dvi 



4nv 
3 

+ 00 



(85) 



When g = g(v), the diffusion tensor (|8ip can be put in the 
form of Eq. (|S"5"]l with 



D \\ = A 



fg 

dv 2 ' 



D X =2A^. 

v dv 



Using Eq. (f55"j) . wc obtain 
8tt . 1 



D H = —A- 



3 v 



+00 



\f(vi)dv 1 + v / vif(vi)dvi 



(86) 



. (87) 



Di_ = —A- 

3 v 



r+oo 

+2v I V\f(vi)dvi 



(88) 
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On the other hand, when h = h(v), the friction term ([52"]) 
can be written as 



I" 1 friction 2Fp / 4^4. V. 

v dv 



(89) 



Using Eq. ([51]). we get 



where 

q(e,t) = \fj a * [2(e-^(r,t))f\ 2 dr, (93) 

is proportional to the phase space volume available to stars 
with an energy less than e. 



1 friction 



2Fp / - 



v f v 

-16nA— f(vi)vldvi. (90) 3.6. Comparison with the two-body encounters theory 

v Jo 



This expression can be obtained directly from Eq. ([521 by 
noting (Binney and Tremaine 2008) that h(v) in Eq. ([53")) 
is similar to the gravitational potential $(r) produced by 
a distribution of mass p(r), where v plays the role of r 
and /(v) the role of p(r). Therefore, if f(v) is isotropic, 
Eq. ([90[) is equivalent to the expression of the gravitational 
field F = -GM(r)r/r 3 produced by a spherically symmet- 
ric distribution of mass, where M(r) is the mass within 
the sphere of radius r. This formula shows that the friction 
is due only to field stars with a velocity less than the ve- 
locity of the test star. This observation was first made by 
Chandrasckhar (1943a). 

The previous expressions for the diffusion and the fric- 
tion coefficients are valid for any isotropic distribution of 
the field particles (of course, when f(v) is the Maxwell dis- 
tribution, we recover the results of Sec. I3.2[) . If we substi- 
tute Eqs. ([57]). ([55]) . and into Eq. g3]), we get a Fokker- 
Planck equation describing the evolution of a test particle in 
a bath with a prescribed distribution /(v)0 Alternatively, 
if we come back to the original Landau kinetic equation 
([28[) . assume an isotropic velocity distribution, and sub- 
stitute the general expressions ([57]) . ([55]) . and ([55]) of the 
diffusion and friction coefficients with now / = f{v,t) we 
obtain the integro-differential equation 



d£ _ 8nA ]_d_ 

dt v 2 dv 



idf (i r 4 



f(vi,t)dvi 



+ SO 



v 1 f(v 1 ,t)dvi ) +/ J f^t^dvi 



(91) 



describing the evolution of the system as a whole. Under 
this form, Eq. (|9ip applies to an artificial infinite homoge- 
neous distribution of stars. This equation has been studied 
by King (1960) in his investigations on the evaporation of 
globular clusters. Within the local approximation, Eq. ([9"T]) 
also represents the simplification of the collision operator 
that occurs in the r.h.s. of the Vlasov-Landau equation ([28]) 
when the velocity distribution of the stars is isotropic. In 
that case, we must restore the space variable and the ad- 
vection term in Eq. (]9ip . From this equation, implement- 
ing an adiabatic approximation, we can derive the orbit- 
averagcd-Fokker-Planck equation which has been used by 
Henon (1961) and Cohn (1980) to study the collisional evo- 
lution of globular clusters. It reads 



dq df_ _ dqdf_ = 8nA ^_ 
de dt dt de de 



J-oo 0€i 



+ 00 



fi de^ 



(92) 



Actually, this description is not self-consistent since the dis- 
tribution f(v) is not steady on the relaxation time in unless it 
is the Maxwell distribution. 



In the previous sections, we have derived the standard ki- 
netic equations of stellar systems from the the Liouville 
equation by using the BBGKY hierarchy. In standard text- 
books of astrophysics (Spitzcr 1987, Binney and Tremaine 
2008), these equations are derived in a different manner. 
One usually starts from the Fokker-Planck equation (|4"3"]l 
and evaluate the diffusion tensor (Av^Av") and the fric- 
tion force (Av 1 - 1 ) by considering the mean effect of a suc- 
cession of two-body encounters. This two-body encounters 
theory was pioneered by Chandrasckhar (1942) and fur- 
ther developed by Rosenbluth et al. (1957). This approach 
directly leads to the expressions ([5T]) and ([82]) of the diffu- 
sion and friction coefficients of a test star in a bath of field 
stars. These expressions are then substituted in the Fokker- 
Planck equation ([43]) . Finally, arguing that the field stars 
and the test star should evolve in the same manner, the 
Fokker-Planck equation is transformed into an integrodif- 
ferential equation ([91]) describing the evolution of the sys- 
tem as a whole (King 1960). When an adiabatic approx- 
imation is implemented (Henon 1961), one finally obtains 
the orbit-averaged- Fokker-Planck equation ([92]) . 

In this paper, we have proceeded the other way round. 
Starting from the Liouville equation, using the BBGKY 
hierarchy, and making a local approximation, we have de- 
rived the Vlasov-Landau equation ([28]) which describes the 
evolution of the system as a whole. Then, making a bath ap- 
proximation^], we have obtained the Fokker-Planck equa- 
tion in the form of Eq. Q36p with the diffusion and friction 
coefficients given by Eqs. ([77]) and ([75[) . Our approach em- 
phasizes the importance of the Landau equation in the ki- 
netic theory of stellar systems, while this equation does not 
appear in the works of Chandrasekhar (1942), Rosenbluth 
et al. (1957), King (1960), and Henon (1961), nor in the 
standard textbooks of stellar dynamics by Spitzer (1987) 
and Binney and Tremaine (2008). 

Actually, the kinetic equation derived by these authors 
is equivalent to the Landau equation but it is written in a 
different form. They write the Fokker-Planck equation in 
the form of Eq. (|4"3"]) with the diffusion coefficient placed 
after the second derivative (d 2 D) while the Landau equa- 
tion ([2"5]) is related to the Fokker-Planck equation ([55]) in 
which the diffusion coefficient is inserted between the first 
derivatives (dDd) . This difference is important on a physi- 
cal point of view for two reasons. First, the Landau equation 
isolates the friction by polarization F po ; while the equa- 
tion derived by Chandrasekhar (1942) and Rosenbluth et 
al. (1957) involves the total friction F friction- Secondly, the 
Landau equation has a nice symmetric structure from which 
we can immediately deduce all the conservation laws of the 
system (conservation of mass, energy, impulse, and angular 
momentum) and the ii-theorem for the Boltzmann entropy 



In the BBGKY hierarchy, this amounts to singling out a 
particular test star and assuming that the distribution of the 
other stars is fixed. 
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(Balescu 2000). These properties are less apparent in the 
equations derived by King (1960) and Henon (1961) for the 
evolution of the system as a whole. It is interesting to note 
that the symmetric structure of the kinetic equation was 
not realized by early stellar dynamicists while the Landau 
equation was known long before in plasma physical 

Finally, the approach based on the Liouville equation 
and on the BBGKY hierarchy is more rigorous than the 
two-body encounters theory because it relaxes the assump- 
tion of locality and does not produce any divergence at large 
scale0. It leads to the generalized Landau equation (j2"4"]l 
that is perfectly well-behaved at large scales contrary to the 
Vlasov-Landau equation (|2"8"|) in which a large-scale cut-off 
has to be introduced in order to avoid a divergence. This di- 
vergence is due to the long-range nature of the gravitational 
potential which precludes a rigorous application of the two- 
body encounters theory that is valid for potentials with 
short-range interactions. Actually, the two-body encounters 
theory is marginally applicable to the gravitational force (it 
only generates a weak logarithmic divergence) and this is 
why it is successful in practice. The generalized Landau 
equation (|24|) represents a conceptual improvement of the 
Vlasov-Landau equation (|2"gj) because it goes beyond the 
local approximation and fully takes into account the spa- 
tial inhomogeneity of the system. Unfortunately, this equa- 
tion is very complicated to be of much practical use. It can 
however be simplified by using angle-action variables as we 
show in the next section. 



4. Kinetic equations with angle-action variables 

4.1. Adiabatic approximation 

In order to deal with spatially inhomogeneous systems, it 
is convenient to introduce angle-action variables (Goldstein 
1956, Binney and Tremaine 2008). Angle-action variables 
have been used by many authors in astrophysics in or- 
der to solve dynamical stability problems (Kalnajs 1977, 
Goodman 1988, Weinberg 1991, Pichon and Cannon 1997, 
Valageas 2006a) or to compute the diffusion and friction co- 
efficients of a test star in a cluster (Lynden-Bell and Kalnajs 
1972, Tremaine and Weinberg 1984, Binney and Lacey 
1988, Weinberg 1998, Nelson and Tremain 1999, Weinberg 
2001, Pichon and Aubert 2006, Valageas 2006b, Chavanis 
2007, Chavanis 2010). By construction, the Hamiltonian H 
in angle and action variables depends only on the actions 
J = ( Ji, Jd) that are constants of the motion; the con- 
jugate coordinates w = (w±, ...,Wd) are called the angles 
(see Appendix [Bj) . Therefore, any distribution of the form 
/ = /(J) is a steady state of the Vlasov equation. According 
to the Jeans theorem, this is not the general form of Vlasov 
steady states. However, if the potential is regular, for all 
practical purposes, any time-independent solution of the 
Vlasov equation may be represented by a distribution of 
the form / = /(J) (strong Jeans theorem). 

We shall assume that the system has reached a quasi sta- 
tionary state (QSS) described by a distribution / = /(J) as 
a result of a violent collisionless relaxation involving only 



23 To our knowledge, the first explicit reference to the Landau 
equation in the astrophysical literature appeared in the paper 
of Kandrup (1981). 

24 It is moreover applicable with almost no modification to 
other systems with long-range interactions for which a two-body 
encounters theory is not justified (Chavanis 2010). 



meanfield effects. Due to finite N effects, the distribution 
function / will slowly evolve in time. Finite N effects are 
taken into account in the "collision" operator appearing in 
the r.h.s. of Eq. (|2"4"|) . Since this term is of order 1/N, the 
effect of "collisions" (granularities, finite N effects, correla- 
tions...) is a very slow process that takes place on a relax- 
ation timescale tn ~ (N/\nN)trj (see Sec. 12. 7[) . Therefore, 
there is a timescale separation between the dynamical time 
to that is the timescale during which the system reaches 
a steady state of the Vlasov equation through phase mix- 
ing and violent collisionless relaxation and the collisional 
relaxation time tn which is the timescale during which the 
system reaches an almost isothermal distribution due to 
finite N effects. 

Because of this timescale separation, the distribution 
function is stationary on the dynamical timescale. It will 
evolve through a sequence of QSSs that are steady states 
of the Vlasov equation, depending only on the actions J, 
slowly changing in time due to the cumulative effect of en- 
counters (finite N effects). Indeed, the system re-adjusts 
itself dynamically at each step of the collisional process. 
The distribution function averaged over a short dynamical 
timescale can be approximated by 



(/(r,v,i))~/(J,t). 



(94) 



Therefore, the distribution function is a function / = 
f(J,t) of the actions only that slowly evolves in time un- 
der the effect of "collisions" . This is similar to an adiabatic 
approximation. The system is approximately in mechanical 
equilibrium at each stage of the dynamics and the "colli- 
sions" slowly drive it towards an almost isothermal distri- 
bution, corresponding to a quasi thermodynamical equilib- 
rium state. 

4.2. Evolution of the system as a whole: a Landau-type 
equation 

Introducing angle-action variables (w,J), the generalized 
Landau equation (|2~4"|) become^ 



r+oo r 

J dr J dwidJ x F(l -> 0)G(t, t - t) 



d£ = d_ 
dt dJ 



F(1^0).A + F( o^i).^_ 



/(J,t)^-(Ji,f), (95) 

TO 



with 

F(l -> 0) = 



(96) 



To obtain Eq. (|9"5)) . we have averaged Eq. (|2"4")l over w (to 
simplify the expressions, the average (.) = (2-7r)~ 3 J dw is 
implicit), written the scalar products as Poisson brackets, 
and used the invariance of the Poisson brackets and of the 



25 The adiabatic assumption (I94|l is consistent with the 
Bogoliubov ansatz of kinetic theory. Since the correlation time of 
the fluctuations is of the order of the dynamical time, or shorter, 
we can freeze the distribution function at time t to compute 
the integral over r (see Section T2.4|l . This distribution function, 
which is a steady state of the Vlasov equation, defines a set 
of angle-action variables (J,w) that we can use to perform the 
integral over r. Then, the distribution function f(J,t) evolves 
with time on a longer timescale according to Eq. (1951) . 
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phase space volume element on a change of canonical vari- It is easy to establish that 
ables. Introducing the Fourier transform of the potential ., /T /tt\ a ^tt\* 

with respect to the angles 



(104) 



A k|kl (J, Ji) = ^ 2d J u[r(J,w)-r(Ji,wi)] 

xe -i(k-w-k 1 -w 1 ) dwdwi) ( Q7 ) 

so that 

U [r(J,w)-r(J 1 ,w 1 )] =^^k,k 1 (J,Ji)e i ( k - w - kl ' Wl \ 

k,ki 



Therefore, the kinetic equation can be rewritten as 



at 



83 



(It 



dwidJi ^2 y^^k,ki(J, Ji) 

k,ki l,h 

xke i(k+l)-w e -i(ki+li)-wi e -i(l-n(J,i)-li-n(Ji,t))r 



f(3,t)^(3 1 ,t). (105) 



we get 
F(l -> 0) 



£A k , kl (J,J 1 )ke i ( k - w - kl - Wl >. 



(98) 
(99) 



Integrating over Wi , and recalling that this expression has 
to be averaged over w, we obtain 

^ = (27r) 3 m 2 ^ • dr J dJrglAk^^Ji)! 2 



k,ki 



^g-sC-k-nOJ.^+ki-ntJi.t^T 



Substituting this expression in Eq. (|95p . wc obtain 



9/ 
8t 



,8_ 
83 



r + co p 

\ dr dwxdJi ^ V] A k ,k! (J, Ji) 











k-^r-ki-^- ) f(3,t)±(3 u t) 



83 



83 1 



f 



(106) 



k,ki l,Ii 



cke < < k,w_kl,Wl )G(t,f-T) 



Al ill (J,Ji)Ie i ( 1 ' w - ll - Wl 



03 



Making the transformation r — » — r, then (k, ki) — > 
(— k, — ki), and adding the resulting expression to Eq. 
([TUB")) , we get 



-i4 Il ,i(Ji,J)lie < ( ll - wl - 1 - w )-^ r 

83 1 



f 



/(J,t)-(Ji,*). 

m 



9J 



-foo 



^ = J(2^) 3 m 2 | f • / dr I d3 t J2 I^WJ, Ji)| s 



k,ki 



(100) 



With angle-action variables, the equations of motion of a 
star determined by the mean field take the very simple form 
(see Appendix [B]) : 

3{t-r) = 3{t) =3, 

w(t - t) = w(t) - fi(J, t)r = w- ft(3, t)r, (101) 



, ke i(k-0(J,t)-k 1 -0(J 1 ,t))T 

f(3,t)^(3 u t). (107) 



v d t 9 

k ■ 83 kl ' aii 

Finally, using the identity (|C.11[) . we obtain the kinetic 
equation 



^ = 7r(2 7 r) 3 mA . J dJik|A kjkl (J, J 



where Jl(J,t) is the angular frequency of the orbit with 
action J. As explained previously, we have neglected the 
variation of the mean field on a timescale of the order of 
the dynamical time so it is considered as "frozen" when we 
compute the stellar trajectories (adiabatic or Bogoliubov 

assumption). Substituting these relations in Eq. ([TOO)) and Tni s kinetic equation was previously derived for systems 



x6[k- Sl(3,t) - ki • fl(3 u t)} 
X l k -^- kl '^) /(J ' <)/(J ^ ) - 



(108) 



making the transformations 1 



-1 and li 



second term (friction term), we obtain successively 
df 8 f + °° f 

— = - m2 ^ ■ dr J dvMJi ^ ^ A kMl (J, J x 

J o J i, i, i i 



k,ki l,li 

x jjgiCk-w-ki-wi) i(l-w(t-T)-li-wi(t-r)) 



A hll {3,3 x )\- 



03 



and 



df 8 r+°° r 

— = -m 2 — •/ dr / dwidJi ^^A k , kl (J,J 

JO J i, u ii. 



k,ki 1,1 

X ] £e i(k-w- ki-wi) e i(l-w-li-wi) e -t(l.n(J,t)-li-n(Ji,t))T 

A hh (3, Jr)l- — - A-u.-^Jx, J)li • — 
a J oJi 

/ 



li in the with arbitrary long-range interactions in various dimen- 
sions of space (Chavanis 2007,2010) and it is here specifi- 
cally applied to stellar systems. Since collective effects are 
neglected, this kinetic equation can be viewed as a Landau- 
type equation with angle-action variables describing the 
evolution of spatially inhomogeneous stellar systems. The 
collisional evolution of these systems is due to a condition of 
resonance k • f2(J, t) = ki • Sl(3i,t) (with (kj., Ji) ^ (k, J)) 
encapsulated in the (5-function. This (5-function expresses 
the conservation of energy. It can be shown (Chavanis 2007) 
j(j ^JL(j 1 f\ (102) that the kinetic equation (|108[) conserves mass M — J f d3 
171 and energy E = J fe(3) d3 and monotonically increases 

the Boltzmann entropy S = — J (/ / to) ln(/ / 'm) d3 (H- 
theorem). However, as explained in Section [2Tfl this equa- 
tion docs not reach a steady state due to the absence of 
statistical equilibrium for stellar systems^. 



A_ h _i(Ji, J)li 



8 
83 ! 



26 For self-gravitating systems in lower dimensions of space, or 
for systems with long-range interactions with a smooth potential 
like the HMF model, a statistical equilibrium state exists. In 
that case, it can be shown that the Landau-type equation (|108|l 
relaxes towards the Boltzmann distribution on a timescale Nto 



xf(3 1 t)^(3 1 ,t). (103) 

provided there are enough resonances (see Chavanis 2007). 
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4.3. Relaxation of a star in a thermal bath: the 
Fokker-Planck equation 

Implementing a test particle approach as in Sec. we find 
that the equation for P(3, t), the probability density of find- 
ing the test star with an action J at time t, is 

n^irfm-^-Y, / <*Jik|A k , kl (J, Jr)| 2 

1_ 1_ *f 



at 



83 



k,ki 

xS\k- «(J) -ki -fi(Ji)] 



(109) 



The angular frequency Q(J) is now a static function deter- 
mined by the distribution /(J) of the field stars. Equation 
(|109[) can be written in the form of a Fokker-Planck equa- 
tion 



dP d ( dP 
— = —— I D^ u —— - PF^ 
dt dJ»\ dJ v po1 

involving a diffusion tensor 

^ Tr(2TT) 3 mJ2 I d3 1 WW\A kM (3,3 1 ) 

x<y[k.n(j)-ki.n(Ji)]/(Ji), 

and a friction by polarization 
F pol = 7r(27r) 3 m ]T [ d3 1 k\A kM (3,3 1 )\ 2 

l, l, . 



(110) 



(111) 



k,k 



Of 



xjpc-ncjj-ki.ncjoik!-^- (Ji). 



(112) 



Writing the Fokker-Planck equation in the usual form 
dP 



d 2 d 
- (D^ U P"\ — ( PF 11 ) 

dt dJW K ' d.L K Motion)' 



with 
, lv _ (AJ»AJ») 



2At 



friction 



(A3) 
At 



(113) 



(114) 



we find that the relation between the friction by polariza- 
tion and the total friction is 



friction pol ' QJ V 



(115) 



Substituting Eqs. ([TTT|) and (jTT2"j) in Eq. ([TT5j) and using an 
integration by parts, we find that the diffusion and friction 
coefficients are given by 

(AJ^AJ V ) 



= 7r(27r) 3 77i I dJi/(Ji)^ WW 

k,ki 

x|A k , kl (j, Ji)| 2 <5(k • n(j) - k x • n(j x )), 



2A< 



(116) 



(AJ) 
At 



7r(27r) 3 m 



/^/(joE^k-^-k.-A 

k.ki 



x|A k)kl (j, Ji)| 2 5(k • n(j) - ki • n(J0). (ii7) 

These expressions can be obtained directly from the 
Hamiltonian equations of motion by expanding the trajec- 
tories of the stars in powers of 1/N in the limit N —> +oo 
(Valageas 2006a). 



Let us assume that the field stars form a thermal bath 
with the Boltzmann distribution 



/(Ji) = Ae- pme(3l \ 



(118) 



where e(J) is the energy of a star in an orbit with action J. 
As we have explained before, this distribution is not defined 
globally for a self-gravitating system. However, it holds ap- 
proximately for stars with low energjQ Using the identity 
de/dJ = il(J) (see Appendix IB"]) . we find that 



dh 
0Ji 



-/3m/(Ji)fi(Ji 



(119) 



Substituting this relation in Eq. (|112p , using the 8- function 
to replace ki ■fl(Ji) by k-f2(J), and comparing the resulting 
expression with Eq. (|111[) , we finally get 

F poi = -D» v {3)pmQ?{3), (120) 

which is the appropriate Einstein relation for our problem. 
For a thermal bath, using Eq. (|120p . the Fokker-Planck 
equation (|110[) can be written as 



dP 
~dt 



d 

tw 1 



D" v {3) 



dP 



■PmPQ?{3\ 



(121) 



where D^(3) is given by Eq. ([TTT1) with Eq. (JTTHJ). 
Recalling that J1(J) = de/83, this equation is similar to 
the Kramers equation in Brownian theory (Kramers 1940). 
This is a drift-diffusion equation describing the evolution 
of the distribution P(3,t) of the test star in an "effective 
potential" U e ff(3) — e(J) produced by the field stars. For 
t — > +oo, the distribution of the test star relaxes towards 
the Boltzmann distribution (|118l) . This takes place on a 
typical relaxation time tn ~ (N/\nN)tu. Again, this is 
valid only in the part of the cluster where the Boltzmann 
distribution holds approximately. 

5. Conclusion 

Starting from the Liouville equation, using a truncation of 
the BBGKY hierarchy at the order 1/N, and neglecting 
collective effects, we have derived a kinetic equation (|2"4"]l in 
physical space that can be viewed as a generalized Landau 
equation. This equation was previously derived by Kandrup 
(1981) using projection operator technics. A nice feature of 
this equation is that it does not present any divergence at 
large scales since the spatial inhomogeneity of the system is 
accounted for. When a local approximation is implemented, 
and a cut-off is introduced heuristically at the Jeans length, 
we recover the Vlasov-Landau equation (|27|) which is the 
standard equation of stellar dynamics. On the other hand, 
using angle-action variables, we have derived a Landau-type 
equation (|108p for fully inhomogeneous stellar systems. We 
have also developed a test particle approach and derived the 
corresponding Fokker-Planck equations f|39[) and (|109[) in 
position-velocity space and angle-action space respectively. 
Explicit expressions have been given for the diffusion and 
friction coefficients. We have distinguished the friction by 
polarization from the total friction. A limitation of the ap- 
proach presented here is that it neglects collective effects. 



27 This description is also valid for self-gravitating systems 
in lower dimension of space, or for other systems with long- 
range interaction, for which a statistical equilibrium state exists 
(Chavanis 2007, 2010). 
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More general kinetic equations, corresponding to Lenard- 
Balescu-type equations taking spatial inhomogeneity and 
collective effects into account, have been derived recently 
by Heyvaerts (2010) from the Liouville equation and by 
Chavanis (2012b) from the Klimontovich equation (these 
approaches based on the BBGKY hierarchy or on the quasi- 
linear approximation are equivalent but the formalism is 
different). These kinetic equations are more general than 
those derived in the present paper, but they are also more 
complicated (to derive and to solve). Therefore, the equa- 
tions presented in this paper may be useful as a first step. 

We have also discussed the differences between the 
present approach based on the BBGKY hierarchy and the 
more classical two-body encounters theory (Chandrasekhar 
1942). The two-body encounters theory, which is usu- 
ally adapted to short-range potentials (Boltzmann 1872, 
Chapman and Cowling 1939), can take strong collisions 
into account so it does not yield any divergence at small 
scales. However, this approach cannot take spatial inhomo- 
geneity into account so it leads to a divergence at large 
scales. This divergence is due to the long-range nature of 
the gravitational potential and the dominance of weak col- 
lisions. In addition, the two-body encounters theory does 
not take collective effects into account; these effects arc spe- 
cific to systems with long-range interactions. By contrast, 
the approach based on the BBGKY hierarchy takes into 
account the spatial inhomogeneity of the system and col- 
lective effects. Therefore, it does not yield any divergence 
at large scales. However, it fails to take strong collisions 
into account due to the weak coupling approximation. In 
a sense, the gravitational potential is intermediate between 
short-range and long-range potentials because both strong 
collisions and weak collisions are relevant. Therefore, the 
approaches adapted to short-range or long-range potentials 
are both marginally applicable (they yield a logarithmic di- 
vergence are large or small scales respectively) . This is why 
the kinetic equations of stellar dynamics can be obtained 
in different manners that turn out to be complementary to 
each other. 

In this paper, we have assumed that the system is iso- 
lated from the surrounding. As a result, the source of noise 
is due to discreteness (finite N) effects internal to the sys- 
tem. The case where the noise is caused by external sources 
(perturbations on a galaxy, cosmological environment on 
dark matter halos...) is also interesting. It has been consid- 
ered by several authors such as Weinberg (2001), Ma and 
Bertschinger (2004), and Pichon and Aubert (2006) who 
developed appropriate kinetic theories. 



Appendix A: The thermodynamic limit 

The kinetic and potential energies in the Hamiltonian 
([1} are comparable provided that v 2 n ~ GNm/R, where 
v m is the root mean square velocity of the stars and R 
the system's size (this scaling can also be obtained from 
the virial theorem). Therefore, the energy scales as E ~ 
Nrnv^ ~ GN 2 m 2 / R and the kinetic temperature, defined 
by k B T = rm4/3, scales as k B T ~ GNm 2 /R ~ E/N. 
The thermodynamic limit of a self-gravitating system cor- 
responds to N — > +oo in such a way that the normalized 
energy e = ERj (GN 2 m 2 ) and the normalized temperature 
r\ = f3GNm 2 / R are of order unity. Of course, the usual 
thermodynamic limit N, V —> +oo with N/V ~ 1 is not 



applicable to self-gravitating systems since these systems 
are spatially inhomogeneous. 

By a suitable normalization of the parameters, we can 
take R ~ 1, m ~ 1, and G ~ 1/N. In this way, E ~ N, 
S ^ N and T ~ 1. This is the proper thermodynamic limit 
for systems with long-range interactions (Kac et al. 1963, 
Campa et al. 2009). We note that the coupling constant G 
scales as 1 /N . The energy and the entropy are extensive but 
they remain fundamentally non-additive. The temperature 
is intensive. The dynamical time tu ~ R/v m ~ 1/y/Gp is 
of order unity (irj ~ 1). 

Other normalizations of the parameters are possible. For 
example, Gilbert (1968) considers the limit TV — s- +oo with 
G ~ 1, R ~ 1, and m ~ 1/N. In that case, E ~ 1, S ~ N, 
T ~ 1/N, and t B ~ 1. On the other hand, de Vega and 
Sanchez (2002) define the thermodynamic limit as N — > 
+oo with m ~ 1, G ~ 1 ; and R ~ N. In that case, E ~ N, 
S ~ N, and T ~ 1. However, the dynamical time tr> ~ 
N diverges with the number of particles. Therefore, this 
normalization may not be convenient to develop the kinetic 
theory of stellar systems. If we impose G ~ 1, E ~ N, T ~ 1 
and t D ~ 1, we get R ~ A 1 / 5 and m ~ N~ 2 / 5 . 

Dimcnsionally, the Jeans length scales as Xj ~ 
{kBT/Gm 2 ^ 1 / 2 , where n is the number density. This is 
the counterpart of the Debye length X B ~ {k B T /ne 2 ) 1 ' 2 
in plasma physics. Since v 2 t ~ GM/R, we find that R ~ 
v m /{Gnra) 1 / 2 ~ [k B T /Gnm 2 ) 1 / 2 ~ Xj. Therefore, the 
Jeans length is of the order of the system's size. 

The gravitational parameter (resp. the plasma parame- 
ter) is defined as the ratio of the interaction strength at the 
mean interparticle distance Gm 2 n x l 3 (resp. e 2 :^ 1 / 3 ) to the 
thermal energy k B T. This leads to g = Gm 2 n 1 / 3 /k B T = 
l/(nA 3 ) 2 / 3 = 1/A 2 / 3 ~ I/Y 2 / 3 (resp. g = e 2 n^ 3 /k B T = 
l/(nX 3 D ) 2 / 3 = I/A 2 / 3 ). Estimating the density by n ~ 
N/R 3 , we find that g ~ rj/N 2 / 3 . Therefore, the expansion 
of the BBGKY hierarchy in terms of the gravitational pa- 
rameter (resp. plasma parameter) g is equivalent to an ex- 
pansion in terms of the inverse of the number of particles 
in the Jeans sphere A = nX 3 ~ N (resp. the inverse of the 
number of particles in the Debye sphere A = nX 3 ^). This 
corresponds to a weak coupling approximation. 

Appendix B: Angle-action variables 

In Section 21 we have explained that, during its colli- 
sional evolution, a stellar system passes by a succession of 
QSSs that are steady states of the Vlasov equation slowly 
changing under the effect of close encounters (finite N ef- 
fects). The slowly varying distribution function /(r,v) de- 
termines a potential <3>(r) and a one-particle Hamiltonian 
e = v 2 /2 + <&(r) that, we assume, is integrable. Therefore, 
it is possible to use angle-action variables constructed with 
this Hamiltonian (Goldstein 1956, Binney and Tremainc 
2008). This construction is done adiabatically, i.e. the dis- 
tribution function, and the angle-action variables, slowly 
change in time. 

A particle with coordinates (r, v) in phase space is de- 
scribed cquivalcntly by the angle-action variables (w, J). 
The Hamiltonian equations for the conjugate variables 
(r, v) are 

dr de dv de _, . . . 
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In terms of the variables (r, v), the dynamics is compli- 
cated because the potential explicitly appears in the sec- 
ond equation. Therefore, this equation dv/dt = — V$ can- 
not be easily integrated except if $ = 0, i.e. for a spa- 
tially homogeneous system. In that case, the velocity v is 
constant and the unperturbed equations of motion reduce 
to r = vt + ro, i.e. to a rectilinear motion at constant 
velocity. Now, the angle-action variables are constructed 
so that the Hamiltonian does not depend on the angles 
w. Therefore, the Hamiltonian equations for the conjugate 
variables (w, J) are 



dw 
~dt 



de 
03 



dJ 

It 



dw 



0. 



(B.2) 



where f2(J) is the angular frequency of the orbit with action 
J. From these equations, we find that J is a constant and 
that w = fi(J)t + wo. Therefore, the equations of motion 
are very simple in these variables. They extend naturally 
the trajectories at constant velocity for spatially homoge- 
neous systems. This is why this choice of variables is rel- 
evant to develop the kinetic theory. Of course, even if the 
description of the motion becomes simple in these variables, 
the complexity of the problem has not completely disap- 
peared. It is now embodied in the relation between position 
and momentum variables and angle and action variables 
which can be quite complicated. 

Appendix C: Calculation of K^" 

Within the local approximation, we can proceed as if the 
system were spatially homogeneous. In that case, the mean 
field force vanishes, (F) = 0, and the unperturbed equa- 
tions of motion (i.e. for N — > +oo) reduce to 



v(i - r) = v(t) = v, 

r(i — r) = r(i) — v(f)r = r vt, 



(C.l) 
(C.2) 



corresponding to a rectilinear motion at constant velocity. 
The collision term in the kinetic equation (|26[) can be writ- 
ten as 



df 
dt 



dv» J \dv» dv\ 



with 



x/(r,v,t)/(r,v 1; t), (C.3) 

K tlv = — f + dr ( dnF^il -> 0,t)F"(l -> 0,f-r). 
m Jo J 

(C.4) 

The force by unit of mass created by particle 1 on particle 
is given by 



du , 

F(1^0) = -m— (r-n), 



(C.5) 



where u(r — r') = — G/\r — r'| is the gravitational potential. 
The Fourier transform, and the inverse Fourier transform, 
of the potential are defined by 



i(k) = J e-*' x u(x) u(x) = J e^u(k) dk. 



(C.6) 



For the gravitational interaction 

(27T) 3 ii(fc) =- — . 



(C.7) 



Substituting Eq. (|C.6I b) in Eq. (|C.5[) . and writing explicitly 
the Lagrangian coordinates, we get 

F(l->0,t-T) = -im J ke lk - {r{t ~ T) - ri{t - T)) u(k) dk. 

(C.8) 

Using the equations of motion (|C.1|) and (|C.2[) . and intro- 
ducing the notations x = r — ri and w = v — Vi , we obtain 



F(l -> 0, t - r) = -im J k e ik ' (x ~ WT) u(k) dk. 



(C.9) 



Therefore 



-m 







dr I rfx / dk I dk'k^kl 



xe l(k+k ' )x e- jk ' WT M(k)M(k'). 
Using the identity 



5(x) 



dk 



(CIO) 
(C.ll) 



and integrating over x and k', we find that 

r+ca /• 

K^ = (2ir) 3 mj dr J dkWk? e lk wr w(fc) 2 . (C.12) 

Performing the transformation r — > — r, then k — > k, and 
adding the resulting expression to Eq. (|C.12[) . we get 

= ^(2ir) 3 m J + dr J dkA^fc*e* -wr fi(fc) a - (C.13) 

Using the identity (jC.lip . we finally obtain 

W = 7r(27r) 3 m J dkk^k w 8{k ■ w)u(fc) 2 , (C.14) 

which leads to Eq. (|27| . 

Introducing a spherical system of coordinates in which 
the z axis is taken in the direction of w, we obtain 



/■+00 

K^ = 7r(27r) 3 m / k 2 
Jo 



2 dk / sin 6d6 



x / d(/)k^k u S(kw cos e)u{kf 



(C.15) 



Using k x = fcsin0cos</>, k y = k sin 8 sine/) and k z = fccos0, 
it is easy to sec that only K xx , K yy and K zz can be non- 
zero. The other components of the matrix K^ v vanish by 
symmetry. Furthermore 



K xx = Kyy = 7r 2 (27r) 3 m / k 2 dk 



x / sind ddk 2 S(kw cosO)u(k) sin 0. 
Jo 

Using the identity S(Xx) = jy-rS(x), we get 



(C.16) 



1 f + 

K xx = Kyy = 7r 2 (27r) 3 m- / k 3 u(k) 2 dk 
w Jo 

x / sin 3 0<5(cos0)d0. (C.17) 
Jo 



20 



Pierre-Henri Chavanis: Kinetic theory of spatially inhomogeneous stellar systems without collective effects 



With the change of variables s = cos 9, we obtain 
1 f + °° 

K xx = K yy = 7r 2 (27r) 3 m— / k 3 u{k) 2 dk 
w Jo 



x / (1 - s 2 )6{s) ds, 



so that, finally, 



K — K 



87r 5 m- / k 3 u(k) 2 dk. 
w Jo 



(C.18) 



(C.19) 



On the other hand, 



1 f^°° 

K zz = 27r 2 (27r) 3 m- / k 3 u(k) 2 dk 
w Jo 

/■7T 

x / sin(9cos 2 6><5(cos(9)d6> = 0. 
Jo 



In conclusion, we obtain 
w 2 5^ v ~ w^w" 



ur 



= A 
with 

r+oo 

A = 87r 5 m / k 3 u(k) 2 dk. 
Jo 



(C.20) 

(C.21) 

(C.22) 



Using Eq. (jCJj) . this leads to Eq. J28 



Appendix D: Another derivation of the Landau 
equation 

For an infinite homogeneous system, the distribution func- 
tion and the two-body correlation function can be written 
as /(0) = /(v,f) and 3(0,1) = g(r — Y\, v, vx, t). In that 
case, Eqs. (fT8f and (fl~9|) become 

% 
dt 



r \ 2 ® f du , . 

(y,t)=m — • / — g(x,v,vx,f)dxdvx, 



(D.l) 



— (X,V,Vx,tJ + w • — (X,V,V lf i) 

du f d d \ „, . f , 
ax \ av avi / m 



(D.2) 



where we have defined x = r — rx and w = v — Vx- Using 
the Bogoliubov ansatz, we shall treat the distribution func- 
tion / as a constant, and determine the asymptotic value 
g(x, v, Vx, +oo) of the correlation function. Introducing the 
Fourier transforms of the potential of interaction and of the 



correlation function, Eq. (jD.ip can be replaced by 
df d f 

-± = (27r) 3 m 2 — • / k{t(fc)lmte(k,v,vx,+oo)l dkdvx, 
dt 9v J 

(D.3) 

where we have used the reality condition g(— k) = <?(k)*. 
On the other hand, taking the Laplace- Fourier transform of 



Eq. (|D.2I) and assuming that no correlation is present ini- 
tially (if there are initial correlations, their effect becomes 
rapidly negligible), we get 

g(k, v, vx , u) = iu{k) — — 

U1{K ■ W — UJ) 

xk -(£-£)« T ' , >s< vi ' , >- (D ' 4) 



Taking the inverse Laplace transform of Eq. (|D.4|) . and us- 
ing the residue theorem, we find that the asymptotic value 
t — > +oo of the correlation function, determined by the pole 
w = 0, is 



g(k, v, vx , +oo) = u(k) 



1 



k • w - i0+ 



xk- (-H--JL)f(v,t)t( Vut ). 
\ov ovi / m 

Using the Plemelj formula 
we get 

Im [g(k, v, vx, +oo)] = 7ru(fc)(5(k ■ w) 



(D.5) 



(D.6) 



(D.7) 



Substituting Eq. (|D.7|) in Eq. (|D.3|) . we obtain the Landau 
equation (|2"7|) . 



Appendix E: Lenard-Balescu equation for 
homogeneous stellar systems 

If we assume that the system is spatially homogeneous (or 
make the local approximation), and take collective effects 
into account, the Vlasov-Landau equation (|2"T|) is replaced 
by the Vlasov-Lenard-Balescu equation 



3]_ df_ N-l^ df 



^+v.^ + ^-(F).^-= 7 r(2 7 r) 3 m-- / W 



d_ 



x5(k ■ w) 



'■(k) 



|e(k,k- v)| 2 V dv" J dv\ 
where e(k, cj) is the dielectric function 



dv 



e(k,w) = 1 + (27r) 3 w(fc) 



k- ^ 
^dv. 



-k- 



(E.2) 



The Landau equation is recovered by taking |e(k, k • v)| 2 = 
1. The Lenard-Balescu equation generalizes the Landau 
equation by replacing the bare potential of interaction u(k) 
by a "dressed" potential of interaction 



Udressedfe, k • v) — 



u{k) 



e(k,k- v) 



(E.3) 



The dielectric function in the denominator takes into ac- 
count the dressing of the particles by their polarization 
cloud. In plasma physics, this term corresponds to a screen- 
ing of the interactions. The Lenard-Balescu equation ac- 
counts for dynamical screening since the velocity v of 
the particles explicitly appears in the effective potential. 
However, for Coulombian interactions, it is a good approxi- 
mation to neglect the deformation of the polarization cloud 
due to the motion of the particles and use the static results 
on screening due to Debye and Hiickel (1923). This amounts 
to replacing the dynamic dielectric function |e(k, k • v)| by 
the static dielectric function |e(k, 0)|. In this approxima- 
tion, Mrf re ssed(k, k ■ v) is replaced by the Debye- Hiickel po- 
tential (2Tr) 3 iiDH(k) = (47re 2 /m 2 )/ \k 2 + k 2 D ) corresponding 
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to Udh{x) = (e 2 /m 2 )e k ° r jr. If we make the same approx- 
imation for stellar systems, we find- tllclt Udressed 

(k, k-v) is 

replaced by 



{2nYu DH {k) = - 



4ttG 



k 2 



(E.4) 



corresponding to ujjh{x) = —Gcos(kjr)/r. In t his a pprox- 
imation, the Vlasov-Lenard-Balescu equation (|E.1[) takes 
the same form as the Vlasov-Landau Eq. (051) except that 
A is now given by A = 2irmG 2 Q with 



l/R 



{k 2 - kj) 2 



dk, 



(E.5) 



where R is the system's size. We see that Q diverges alge- 
braically, as (A,/ — R)~ l . when R — > Xj instead of yield- 
ing a finite Coulomb factor In A ~ \nN when collective 
effects are neglected 2 ^!. This naive approach shows that col- 
lective effects tend to increase the diffusion coefficient or 
tend to reduce the relaxation time. This is the conclusion 
reached by Weinberg (1993) with a more precise approach. 
However, this approach is not fully satisfactory since the 
system is assumed to be spatially homogeneous and the or- 
dinary Lenard-Balescu equation is used. In that case, the 
divergence when R — > Xj is a manifestation of the Jeans in- 
stability that a spatially homogeneous self-gravitating sys- 
tem experiences when the size of the perturbation over- 
comes the Jeans length. For inhomogeneous systems, the 
Jeans instability is suppressed so the results of Weinberg 
should be used with caution. They suggest, however, an 
increase of the relaxation time for an inhomogeneous stel- 
lar system that is not far from being unstable. Hcyvaerts 
(2010) derived a more satisfactory Lenard-Balescu equation 
that is valid for spatially inhomogeneous self-gravitating 
systems. This equation does not present any divergence at 
large scales. However, this equation is complicated and it 
is difficult to measure the importance of collective effects. 

Appendix F: Multi-species systems 

It is straightforward to generalize the kinetic theory of stel- 
lar systems for several species. The Vlasov-Landau equation 
is replaced by 



df a ,„ ^ d 
dt v ' dvf 



^n(k.w)!i 2 (fc) 



(F.l) 



where / a (r,v,t) is the distribution function of species a 
normalized such that J f a drdv ~ N a m a . We can use this 
equation to give a new interpretation of the test particle 
approach developed in Sec. |31 We make three assumptions: 
(i) We assume that the system is composed of two types 
of stars, the test stars with mass m and the field stars 
with mass To/; (ii) we assume that the number of test stars 
is much lower than the number of field stars; (iii) we as- 
sume that the field stars are in a steady distribution /(r, v). 



In plasma 



physics, on the contrary, Q = 
f* L k s /[(k 2 + k 2 D ) 2 ]dk is well-behaved as k -> 0. When 
collective effects are taken into account, there is no divergence 
at large scales and the Debye length appears naturally. In the 
dominant approximation Q ~ In A with A = n\ A D . 



Because of assumption (ii), the collisions between field stars 
and test stars are negligible so that the field stars remain 
in their steady state. The collisions between test stars are 
also negligible, so they only evolve due to collisions with 
the field stars. Therefore, if we call P(r, v, t) the distribu- 
tion function of the test stars (to have notations similar to 
those of Sec. [3] with, however, a different interpretation), 
its evolution is given by the Fokker-Planck equation 

dP 8 f 

— = n{2nf-^ J k»k»6(k-w)u 2 (k) 

X ( m//l ^" WP ^l) dVldk - (K2) 
The diffusion and friction coefficients are given by 

D*" = n{2nfm f J k>*k"5{k ■ w)w(/c) 2 /i dvidk, (F.3) 



F pol = A^fm f fc"r<5(k • w)u(fc) 2 |4 dv i rfk - ( F 



4) 



We note that the diffusion is due to the fluctuations of the 
gravitational force produced by the field stars, while the 
friction by polarization is due to the perturbation on the 
distribution of the field stars caused by the test stars. This 
explains the occurrence of the masses to/ and m in Eqs. 
(|F.3[) and (|F.4j) respectively. 

Using Eq. (|45|) and noting that 



we get 



friction 



1 



TO/ 

m 



pol • 



(F.5) 



(F.6) 



If we assume furthermore that to 3> to/, we find that 
F friction — ^ pol ■ However, in general, the friction force is 
different from the friction by polarization. The other results 
of Sec.[3]can be easily generalized to multi-species systems. 



Appendix G: Temporal correlation tensor of the 
gravitational force 

The diffusion of the stars is due to the fluctuations of 
the gravitational force. For an infinite homogeneous sys- 
tem (or in the local approximation), the diffusion tensor 
can be derived from the formula (|42[) that is well-known 
in Brownian theory. This expression involves the tempo- 
ral auto-correlation tensor of the gravitational force expe- 
rienced by a star. It can be written as 

{F^{t)F"{t - t)) = — [ dndv! -> 0, t) 

TO J 

xHHO.i-rj/W. (G.l) 
Let us first compute the tensor 

Q» v = — I -> 0,^(1 -> 0,t-T)dr v (G.2) 

TO J 

Proceeding as in Appendix [Cj we find 



Q» v = (2tt) 3 to J dki^Ar v e <k - WT u(jfc) a 



(G.3) 
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Introducing a system of spherical coordinates with the z- 
axis in the direction of w, and using Eq. (|C.7[) . we obtain 
after some calculations 



where / is the three-dimensional Fourier transform of the 
distribution function. For the Maxwell distribution (j50|) . we 
obtain 



= 2nmG' 



1 w 2 8^ u - w^w" 



(G.4) (F»{t)F»{t-r)) = (2irfmp 



According to Eq. (jC.4[) , we have 
Q"" dr. 



(G.ll) 



(G.5) 



Eqs. (|G.4[) and (|G.5|I leads to Eq. (J29J) with In A replaced 
by 



InA' =/*-*:. 



(G.6) 



where t m i n and t max are appropriate cut-offs. The upper 
cut-off should be identified with the dynamical time t d ■ On 
the other hand, the divergence at short times is due to the 
inadequacy of our assumption of straight-line trajectories 
to describe very close encounters. If we take i m j„ = Ai,/u m 
and t max = Xj/v m , we find that In A' = In A. In Appendix 
[Cl we have calculated K^ v by integrating first over time 
then over space. This yields the logarithmic factor (l30l) . 
Here, we have integrated first over space then over time. 
This yields the logarithmic factor (|G.6j) . As discussed by 
Lee (1968), these two approaches are essentially equivalent. 
We remark, however, that the calculations of Sec. [C] can be 
performed for arbitrary potentials, while the calculations of 
this section explicitly use the specific form of the gravita- 
tional potential. 

According to Eqs 



Introducing a system of spherical coordinates with the z- 
axis in the direction of v, and using Eq. (|C.7|) . we obtain 
after some calculations 



(G.12) 



where G^ v (x) is defined in Section [5751 In particular, 



(G.13) 



Integrating Eq. (|G.12[) over time, we recover the expression 
(|6"5)) of the diffusion tensor for a Maxwellian distribution 
with In A' instead of In A. 

Finally, the auto-correlation tensor of the gravitational 
field at two different points (at the same time) is 



(F M (0)F"(r)) = 2irnm 2 G 



(G.14) 



(tGTl) . (tCL2]l and <l(L4l) 

auto-correlation function can be written as 



the force References 



(F"(t)F v (t - t)) = 2irmG 2 ^- J dvi 



S^w 2 - w^w" 



/(vi) 



(G.7) 



In particular 

(F(t) ■ F(t - t)) = 47rmG 2 - [ , /(Vl) , dvi- (G.8) 

tJ |v-vi| 

The t _1 decay of the auto-correlation function of the grav- 
itational force was first derived by Chandrasekhar (1944) 
with a different method. This result has also been obtained, 
and discussed, by Cohen et al. (1950) and Lee (1968). 
According to Eqs. (l4"2")l and (|G.7|) . the diffusion tensor is 
given bjo 



= 27rmG 2 lnA' / dw 1 



S^w 2 



/(vi) 



(G.9) 



This returns Eq. (|77|) with In A' instead of In A. 

When /(vi) is the Maxwell distribution (|5"0"|) . we can 
compute the force auto-correlation tensor from Eq. (|G.7[) 
by using the Rosenbluth potentials as in Section 13.51 
Alternatively, combining Eqs. (|G.1[) and (|G.3[) . we have 

{F»(t)F v (t-T)) =(2tt) 6 to J dk^fc !y e lk - vr {t(fc) 2 /(kr), 

(G.10) 



Actually, when the force auto-correlation function decreases 
as t~ , Eq. (|42p is not valid anymore and ((Av) 2 ) behaves as 
tint (Lee 1968). 
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